pkg/seaice/seaice_solve4temp.F

9637aec598 Jean*0001 #include "SEAICE_OPTIONS.h"
1043a55cc1 Jean*0002 #ifdef ALLOW_EXF
                0003 # include "EXF_OPTIONS.h"
                0004 #endif
772b2ed80e Gael*0005 #ifdef ALLOW_AUTODIFF
                0006 # include "AUTODIFF_OPTIONS.h"
                0007 #endif
9637aec598 Jean*0008
                0009 CBOP
                0010 C     !ROUTINE: SEAICE_SOLVE4TEMP
                0011 C     !INTERFACE:
                0012       SUBROUTINE SEAICE_SOLVE4TEMP(
                0013      I   UG, HICE_ACTUAL, HSNOW_ACTUAL,
840c7fba30 Gael*0014 #ifdef SEAICE_CAP_SUBLIM
52ff14d141 Ian *0015      I   F_lh_max,
                0016 #endif
f5282c5b03 Gael*0017      I   TSURFin,
                0018      O   TSURFout,
4dd39c50d9 Mart*0019 #ifdef SEAICE_USE_GROWTH_ADX
                0020      O   F_io_net, F_ia_net,
                0021 #endif /* SEAICE_USE_GROWTH_ADX */
a676451ac2 Jean*0022      O   F_ia, IcePenetSW,
69d8fa39b5 Mart*0023      O   FWsublim,
9637aec598 Jean*0024      I   bi, bj, myTime, myIter, myThid )
                0025
                0026 C     !DESCRIPTION: \bv
                0027 C     *==========================================================*
                0028 C     | SUBROUTINE SOLVE4TEMP
                0029 C     | o Calculate ice growth rate, surface fluxes and
                0030 C     |   temperature of ice surface.
                0031 C     |   see Hibler, MWR, 108, 1943-1973, 1980
                0032 C     *==========================================================*
                0033 C     \ev
                0034
                0035 C     !USES:
                0036       IMPLICIT NONE
                0037 C     === Global variables ===
                0038 #include "SIZE.h"
                0039 #include "GRID.h"
                0040 #include "EEPARAMS.h"
c68247b0e4 Jean*0041 #include "PARAMS.h"
9637aec598 Jean*0042 #include "FFIELDS.h"
ccaa3c61f4 Patr*0043 #include "SEAICE_SIZE.h"
9637aec598 Jean*0044 #include "SEAICE_PARAMS.h"
ccaa3c61f4 Patr*0045 #include "SEAICE.h"
9637aec598 Jean*0046 #include "DYNVARS.h"
                0047 #ifdef ALLOW_EXF
                0048 # include "EXF_FIELDS.h"
                0049 #endif
6d55e2b45a Mart*0050 #ifdef ALLOW_AUTODIFF_TAMC
                0051 # include "tamc.h"
                0052 #endif
9637aec598 Jean*0053
a676451ac2 Jean*0054 C     !INPUT PARAMETERS:
                0055 C     UG           :: atmospheric wind speed (m/s)
9637aec598 Jean*0056 C     HICE_ACTUAL  :: actual ice thickness
                0057 C     HSNOW_ACTUAL :: actual snow thickness
4dd39c50d9 Mart*0058 C     TSURFin    :: surface temperature of ice/snow in Kelvin
a676451ac2 Jean*0059 C     bi,bj      :: tile indices
                0060 C     myTime     :: current time in simulation
                0061 C     myIter     :: iteration number in simulation
                0062 C     myThid     :: my Thread Id number
                0063 C     !OUTPUT PARAMETERS:
4dd39c50d9 Mart*0064 C     TSURFout   :: updated surface temperature of ice/snow in Kelvin
                0065 C     F_io_net   :: upward conductive heat flux through seaice+snow
                0066 C     F_ia_net   :: net heat flux divergence at the sea ice/snow surface:
                0067 C                 includes ice conductive fluxes and atmospheric fluxes (W/m^2)
a676451ac2 Jean*0068 C     F_ia       :: upward seaice/snow surface heat flux to atmosphere (W/m^2)
                0069 C     IcePenetSW :: short wave heat flux transmitted through ice (+=upward)
                0070 C     FWsublim   :: fresh water (mass) flux due to sublimation (+=up)(kg/m^2/s)
6e70381b89 Jean*0071 C---- Notes:
                0072 C     1) should add IcePenetSW to F_ia to get the net surface heat flux
                0073 C        from the atmosphere (IcePenetSW not currently included in F_ia)
                0074 C     2) since zero ice/snow heat capacity is assumed, all the absorbed Short
                0075 C        -Wave is used to warm the ice/snow surface (heating profile ignored).
                0076 C----------
a676451ac2 Jean*0077       _RL UG          (1:sNx,1:sNy)
                0078       _RL HICE_ACTUAL (1:sNx,1:sNy)
                0079       _RL HSNOW_ACTUAL(1:sNx,1:sNy)
840c7fba30 Gael*0080 #ifdef SEAICE_CAP_SUBLIM
a676451ac2 Jean*0081       _RL F_lh_max    (1:sNx,1:sNy)
52ff14d141 Ian *0082 #endif
f5282c5b03 Gael*0083       _RL TSURFin     (1:sNx,1:sNy)
                0084       _RL TSURFout    (1:sNx,1:sNy)
4dd39c50d9 Mart*0085       _RL F_io_net    (1:sNx,1:sNy)
                0086       _RL F_ia_net    (1:sNx,1:sNy)
a676451ac2 Jean*0087       _RL F_ia        (1:sNx,1:sNy)
                0088       _RL IcePenetSW  (1:sNx,1:sNy)
                0089       _RL FWsublim    (1:sNx,1:sNy)
9637aec598 Jean*0090       INTEGER bi, bj
                0091       _RL     myTime
                0092       INTEGER myIter, myThid
1043a55cc1 Jean*0093 CEOP
9637aec598 Jean*0094
1043a55cc1 Jean*0095 #if defined(ALLOW_ATM_TEMP) && defined(ALLOW_DOWNWARD_RADIATION)
9637aec598 Jean*0096 C     !LOCAL VARIABLES:
                0097 C     === Local variables ===
c68247b0e4 Jean*0098 C     i, j  :: Loop counters
840c7fba30 Gael*0099 C     kSurface  :: vertical index of surface layer
9637aec598 Jean*0100       INTEGER i, j
840c7fba30 Gael*0101       INTEGER kSurface
9637aec598 Jean*0102       INTEGER ITER
999f896fe3 Mart*0103       _RL D1, D1I
                0104       _RL D3(1:sNx,1:sNy)
9b4e2be4e7 Mart*0105       _RL TMELT, XKI, XKS, HCUT, recip_HCUT, XIO
b58e51ce4e Jean*0106 C     SurfMeltTemp :: Temp (K) above which wet-albedo values are used
2192eec603 Mart*0107       _RL SurfMeltTemp
a676451ac2 Jean*0108 C     effConduct :: effective conductivity of combined ice and snow
2192eec603 Mart*0109       _RL effConduct(1:sNx,1:sNy)
a676451ac2 Jean*0110 C     lhSublim :: latent heat of sublimation (SEAICE_lhEvap + SEAICE_lhFusion)
                0111       _RL lhSublim
                0112 C     t1,t2,t3,t4 :: powers of temperature
                0113       _RL  t1, t2, t3, t4
9637aec598 Jean*0114
6e70381b89 Jean*0115 C-    Constants to calculate Saturation Vapor Pressure
aa0478ba9b Jean*0116 C     Maykut Polynomial Coeff. for Sat. Vapor Press
a676451ac2 Jean*0117       _RL C1, C2, C3, C4, C5, QS1
6e70381b89 Jean*0118 C     Extended temp-range expon. relation Coeff. for Sat. Vapor Press
a676451ac2 Jean*0119       _RL lnTEN
9637aec598 Jean*0120       _RL aa1,aa2,bb1,bb2,Ppascals,cc0,cc1,cc2,cc3t
                0121 C     specific humidity at ice surface variables
a676451ac2 Jean*0122       _RL mm_pi,mm_log10pi
9637aec598 Jean*0123
4dd39c50d9 Mart*0124 C     tempFrz  :: ocean temperature in contact with ice
                0125 C                 (=seawater freezing point) (K)
be38bd5ad7 Jean*0126 C     F_c      :: conductive heat flux through seaice+snow (+=upward)
aa0478ba9b Jean*0127 C     F_lwu    :: upward long-wave surface heat flux (+=upward)
                0128 C     F_sens   :: sensible surface heat flux         (+=upward)
a676451ac2 Jean*0129 C     F_lh     :: latent heat flux (sublimation) (+=upward)
6e70381b89 Jean*0130 C     qhice    :: saturation vapor pressure of snow/ice surface
                0131 C     dqh_dTs  :: derivative of qhice w.r.t snow/ice surf. temp
                0132 C     dFia_dTs :: derivative of surf heat flux (F_ia) w.r.t surf. temp
4dd39c50d9 Mart*0133       _RL tempFrz    (1:sNx,1:sNy)
aa0478ba9b Jean*0134       _RL F_c        (1:sNx,1:sNy)
a676451ac2 Jean*0135       _RL F_lwu      (1:sNx,1:sNy)
                0136       _RL F_sens     (1:sNx,1:sNy)
                0137       _RL F_lh       (1:sNx,1:sNy)
                0138       _RL qhice      (1:sNx,1:sNy)
6e70381b89 Jean*0139       _RL dqh_dTs    (1:sNx,1:sNy)
aa0478ba9b Jean*0140       _RL dFia_dTs   (1:sNx,1:sNy)
a676451ac2 Jean*0141       _RL absorbedSW (1:sNx,1:sNy)
                0142       _RL penetSWFrac
112191cc93 Jean*0143       _RL delTsurf
136908bfac Ian *0144
a676451ac2 Jean*0145 C     local copies of global variables
                0146       _RL tsurfLoc   (1:sNx,1:sNy)
112191cc93 Jean*0147       _RL tsurfPrev  (1:sNx,1:sNy)
a676451ac2 Jean*0148       _RL atempLoc   (1:sNx,1:sNy)
                0149       _RL lwdownLoc  (1:sNx,1:sNy)
                0150       _RL ALB        (1:sNx,1:sNy)
                0151       _RL ALB_ICE    (1:sNx,1:sNy)
                0152       _RL ALB_SNOW   (1:sNx,1:sNy)
                0153 C     iceOrNot :: this is HICE_ACTUAL.GT.0.
                0154       LOGICAL iceOrNot(1:sNx,1:sNy)
7c50f07931 Mart*0155 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0156 C     ikey :: tape key (depends on iteration)
                0157       INTEGER ikey
7c50f07931 Mart*0158 #endif
a676451ac2 Jean*0159 #ifdef SEAICE_DEBUG
                0160 #endif /* SEAICE_DEBUG */
                0161
                0162 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
9637aec598 Jean*0163
6d55e2b45a Mart*0164 #ifdef ALLOW_AUTODIFF_TAMC
1574069d50 Mart*0165 C     Since this is a "local" tape, is does not need to include number
                0166 C     of threads or tiles per mpi-process. It does not even lead to a
                0167 C     common block.
edb6656069 Mart*0168 CADJ INIT loctape_solve4temp = COMMON, NMAX_TICE
6d55e2b45a Mart*0169 #endif /* ALLOW_AUTODIFF_TAMC */
                0170
6e70381b89 Jean*0171 C-    MAYKUT CONSTANTS FOR SAT. VAP. PRESSURE TEMP. POLYNOMIAL
2192eec603 Mart*0172       C1=    2.7798202  _d -06
                0173       C2=   -2.6913393  _d -03
                0174       C3=    0.97920849 _d +00
                0175       C4= -158.63779    _d +00
                0176       C5= 9653.1925     _d +00
9637aec598 Jean*0177       QS1=0.622 _d +00/1013.0 _d +00
6e70381b89 Jean*0178 C-    Extended temp-range expon. relation Coeff. for Sat. Vapor Press
a676451ac2 Jean*0179       lnTEN = LOG(10.0 _d 0)
9637aec598 Jean*0180       aa1 = 2663.5 _d 0
                0181       aa2 = 12.537 _d 0
                0182       bb1 = 0.622 _d 0
2192eec603 Mart*0183       bb2 = 1.0 _d 0 - bb1
9637aec598 Jean*0184       Ppascals = 100000. _d 0
2192eec603 Mart*0185 C     cc0 = TEN ** aa2
a676451ac2 Jean*0186       cc0 = EXP(aa2*lnTEN)
2192eec603 Mart*0187       cc1 = cc0*aa1*bb1*Ppascals*lnTEN
9637aec598 Jean*0188       cc2 = cc0*bb2
                0189
0320e25227 Mart*0190       IF ( usingPCoords ) THEN
840c7fba30 Gael*0191        kSurface        = Nr
                0192       ELSE
                0193        kSurface        = 1
                0194       ENDIF
9637aec598 Jean*0195
                0196 C     SENSIBLE HEAT CONSTANT
fff6be1885 Mart*0197       D1=SEAICE_dalton*SEAICE_cpAir*SEAICE_rhoAir
9637aec598 Jean*0198
                0199 C     ICE LATENT HEAT CONSTANT
136908bfac Ian *0200       lhSublim = SEAICE_lhEvap + SEAICE_lhFusion
                0201       D1I=SEAICE_dalton*lhSublim*SEAICE_rhoAir
9637aec598 Jean*0202
                0203 C     MELTING TEMPERATURE OF ICE
2192eec603 Mart*0204       TMELT        = celsius2K
9637aec598 Jean*0205
                0206 C     ICE CONDUCTIVITY
                0207       XKI=SEAICE_iceConduct
                0208
                0209 C     SNOW CONDUCTIVITY
                0210       XKS=SEAICE_snowConduct
                0211
                0212 C     CUTOFF SNOW THICKNESS
a676451ac2 Jean*0213 C     Snow-Thickness above HCUT: SW optically thick snow (=> snow-albedo).
                0214 C     Snow-Thickness below HCUT: linear transition to ice-albedo
                0215       HCUT = SEAICE_snowThick
28e63a2c6a Jean*0216       recip_HCUT = 0. _d 0
                0217       IF ( HCUT.GT.0. _d 0 ) recip_HCUT = 1. _d 0 / HCUT
9637aec598 Jean*0218
                0219 C     PENETRATION SHORTWAVE RADIATION FACTOR
                0220       XIO=SEAICE_shortwave
                0221
b58e51ce4e Jean*0222 C     Temperature Threshold for wet-albedo:
                0223       SurfMeltTemp = TMELT + SEAICE_wetAlbTemp
6e70381b89 Jean*0224 C     old SOLVE4TEMP_LEGACY setting, consistent with former celsius2K value:
                0225 c     TMELT        = 273.16  _d +00
                0226 c     SurfMeltTemp = 273.159 _d +00
                0227
c68247b0e4 Jean*0228 C     Initialize variables
9637aec598 Jean*0229       DO J=1,sNy
2192eec603 Mart*0230        DO I=1,sNx
873b79f1c8 Jean*0231 C     initialise output arrays:
                0232         TSURFout (I,J) = TSURFin(I,J)
                0233         F_ia     (I,J) = 0. _d 0
4dd39c50d9 Mart*0234         F_ia_net (I,J) = 0. _d 0
                0235         F_io_net (I,J) = 0. _d 0
873b79f1c8 Jean*0236         IcePenetSW(I,J)= 0. _d 0
                0237         FWsublim (I,J) = 0. _d 0
2192eec603 Mart*0238 C     HICE_ACTUAL is modified in this routine, but at the same time
                0239 C     used to decided where there is ice, therefore we save this information
                0240 C     here in a separate array
a676451ac2 Jean*0241         iceOrNot  (I,J) = HICE_ACTUAL(I,J) .GT. 0. _d 0
                0242         absorbedSW(I,J) = 0. _d 0
2192eec603 Mart*0243         qhice    (I,J) = 0. _d 0
6e70381b89 Jean*0244         dqh_dTs  (I,J) = 0. _d 0
69d8fa39b5 Mart*0245         F_lh     (I,J) = 0. _d 0
2192eec603 Mart*0246         F_lwu    (I,J) = 0. _d 0
                0247         F_sens   (I,J) = 0. _d 0
aa0478ba9b Jean*0248 C     Make a local copy of LW, surface & atmospheric temperatures
f5282c5b03 Gael*0249         tsurfLoc (I,J) = TSURFin(I,J)
                0250 c       tsurfLoc (I,J) = MIN( celsius2K+MAX_TICE, TSURFin(I,J) )
aa0478ba9b Jean*0251         lwdownLoc(I,J) = MAX( MIN_LWDOWN, LWDOWN(I,J,bi,bj) )
6e70381b89 Jean*0252         atempLoc (I,J) = MAX( celsius2K+MIN_ATEMP, ATEMP(I,J,bi,bj) )
9637aec598 Jean*0253
840c7fba30 Gael*0254 c     FREEZING TEMP. OF SEA WATER (K)
                0255         tempFrz(I,J) = SEAICE_dTempFrz_dS *salt(I,J,kSurface,bi,bj)
                0256      &     + SEAICE_tempFrz0 + celsius2K
28e63a2c6a Jean*0257
aa0478ba9b Jean*0258 C     Now determine fixed (relative to tsurf) forcing term in heat budget
                0259
999f896fe3 Mart*0260         IF(HSNOW_ACTUAL(I,J).GT.0.0) THEN
a676451ac2 Jean*0261 C     Stefan-Boltzmann constant times emissivity
999f896fe3 Mart*0262          D3(I,J)=SEAICE_snow_emiss*SEAICE_boltzmann
                0263 #ifdef EXF_LWDOWN_WITH_EMISSIVITY
                0264 C     This is now [(1-emiss)*lwdown - lwdown]
aa0478ba9b Jean*0265          lwdownLoc(I,J) = SEAICE_snow_emiss*lwdownLoc(I,J)
999f896fe3 Mart*0266 #else /* use the old hard wired inconsistent value  */
aa0478ba9b Jean*0267          lwdownLoc(I,J) = 0.97 _d 0*lwdownLoc(I,J)
999f896fe3 Mart*0268 #endif /* EXF_LWDOWN_WITH_EMISSIVITY */
                0269         ELSE
a676451ac2 Jean*0270 C     Stefan-Boltzmann constant times emissivity
999f896fe3 Mart*0271          D3(I,J)=SEAICE_ice_emiss*SEAICE_boltzmann
                0272 #ifdef EXF_LWDOWN_WITH_EMISSIVITY
                0273 C     This is now [(1-emiss)*lwdown - lwdown]
aa0478ba9b Jean*0274          lwdownLoc(I,J) = SEAICE_ice_emiss*lwdownLoc(I,J)
999f896fe3 Mart*0275 #else /* use the old hard wired inconsistent value  */
aa0478ba9b Jean*0276          lwdownLoc(I,J) = 0.97 _d 0*lwdownLoc(I,J)
999f896fe3 Mart*0277 #endif /* EXF_LWDOWN_WITH_EMISSIVITY */
                0278         ENDIF
2192eec603 Mart*0279        ENDDO
9637aec598 Jean*0280       ENDDO
                0281
                0282       DO J=1,sNy
2192eec603 Mart*0283        DO I=1,sNx
9637aec598 Jean*0284
                0285 C     DECIDE ON ALBEDO
2192eec603 Mart*0286         IF ( iceOrNot(I,J) ) THEN
4fecc0d4b4 Jean*0287
2192eec603 Mart*0288          IF ( YC(I,J,bi,bj) .LT. 0.0 _d 0 ) THEN
                0289           IF (tsurfLoc(I,J) .GE. SurfMeltTemp) THEN
                0290            ALB_ICE (I,J)   = SEAICE_wetIceAlb_south
                0291            ALB_SNOW(I,J)   = SEAICE_wetSnowAlb_south
                0292           ELSE                  ! no surface melting
                0293            ALB_ICE (I,J)   = SEAICE_dryIceAlb_south
                0294            ALB_SNOW(I,J)   = SEAICE_drySnowAlb_south
                0295           ENDIF
                0296          ELSE                   !/ Northern Hemisphere
                0297           IF (tsurfLoc(I,J) .GE. SurfMeltTemp) THEN
                0298            ALB_ICE (I,J)   = SEAICE_wetIceAlb
                0299            ALB_SNOW(I,J)   = SEAICE_wetSnowAlb
                0300           ELSE                  ! no surface melting
                0301            ALB_ICE (I,J)   = SEAICE_dryIceAlb
                0302            ALB_SNOW(I,J)   = SEAICE_drySnowAlb
                0303           ENDIF
                0304          ENDIF                  !/ Albedo for snow and ice
9637aec598 Jean*0305
a676451ac2 Jean*0306 C     If actual snow thickness exceeds the cutoff thickness, use snow albedo
2192eec603 Mart*0307          IF (HSNOW_ACTUAL(I,J) .GT. HCUT) THEN
                0308           ALB(I,J) = ALB_SNOW(I,J)
a676451ac2 Jean*0309          ELSEIF ( HCUT.LE.ZERO ) THEN
                0310           ALB(I,J) = ALB_ICE(I,J)
2192eec603 Mart*0311          ELSE
a676451ac2 Jean*0312 C     otherwise, use linear transition between ice and snow albedo
9b4e2be4e7 Mart*0313           ALB(I,J) = MIN( ALB_ICE(I,J) + HSNOW_ACTUAL(I,J)*recip_HCUT
a676451ac2 Jean*0314      &                                 *(ALB_SNOW(I,J) -ALB_ICE(I,J))
                0315      &                  , ALB_SNOW(I,J) )
2192eec603 Mart*0316          ENDIF
9637aec598 Jean*0317
a676451ac2 Jean*0318 C     Determine the fraction of shortwave radiative flux remaining
                0319 C     at ocean interface after scattering through the snow and ice.
                0320 C     If snow is present, no radiation penetrates through snow+ice
2192eec603 Mart*0321          IF (HSNOW_ACTUAL(I,J) .GT. 0.0 _d 0) THEN
a676451ac2 Jean*0322           penetSWFrac = 0.0 _d 0
2192eec603 Mart*0323          ELSE
a676451ac2 Jean*0324           penetSWFrac = XIO*EXP(-1.5 _d 0 * HICE_ACTUAL(I,J))
2192eec603 Mart*0325          ENDIF
a676451ac2 Jean*0326 C     The shortwave radiative flux leaving ocean beneath ice (+=up).
                0327          IcePenetSW(I,J) = -(1.0 _d 0 - ALB(I,J))
                0328      &                    *penetSWFrac * SWDOWN(I,J,bi,bj)
                0329 C     The shortwave radiative flux convergence in the seaice.
                0330          absorbedSW(I,J) =  (1.0 _d 0 - ALB(I,J))
                0331      &        *(1.0 _d 0 - penetSWFrac)* SWDOWN(I,J,bi,bj)
9637aec598 Jean*0332
6e70381b89 Jean*0333 C     The effective conductivity of the two-layer snow/ice system.
aa0478ba9b Jean*0334 C     Set a minimum sea ice thickness of 5 cm to bound
a676451ac2 Jean*0335 C     the magnitude of conductive heat fluxes.
                0336 Cif   * now taken care of by SEAICE_hice_reg in seaice_growth
                0337 c        hice_tmp = max(HICE_ACTUAL(I,J),5. _d -2)
2192eec603 Mart*0338          effConduct(I,J) = XKI * XKS /
a676451ac2 Jean*0339      &        (XKS * HICE_ACTUAL(I,J) + XKI * HSNOW_ACTUAL(I,J))
9637aec598 Jean*0340
                0341 #ifdef SEAICE_DEBUG
a676451ac2 Jean*0342          IF ( (I .EQ. SEAICE_debugPointI) .AND.
52ff14d141 Ian *0343      &        (J .EQ. SEAICE_debugPointJ) ) THEN
2192eec603 Mart*0344           print '(A,i6)','-----------------------------------'
                0345           print '(A,i6)','ibi merged initialization ', myIter
                0346           print '(A,i6,4(1x,D24.15))',
                0347      &         'ibi iter, TSL, TS     ',myIter,
f5282c5b03 Gael*0348      &         tsurfLoc(I,J), TSURFin(I,J)
2192eec603 Mart*0349           print '(A,i6,4(1x,D24.15))',
                0350      &         'ibi iter, TMELT       ',myIter,TMELT
                0351           print '(A,i6,4(1x,D24.15))',
                0352      &         'ibi iter, HIA, EFKCON ',myIter,
                0353      &         HICE_ACTUAL(I,J), effConduct(I,J)
                0354           print '(A,i6,4(1x,D24.15))',
                0355      &         'ibi iter, HSNOW       ',myIter,
                0356      &         HSNOW_ACTUAL(I,J), ALB(I,J)
                0357           print '(A,i6)','-----------------------------------'
                0358           print '(A,i6)','ibi energy balance iterat ', myIter
                0359          ENDIF
300dea8deb Jean*0360 #endif /* SEAICE_DEBUG */
4fecc0d4b4 Jean*0361
2192eec603 Mart*0362         ENDIF                   !/* iceOrNot */
                0363        ENDDO                    !/* i */
                0364       ENDDO                     !/* j */
a676451ac2 Jean*0365
                0366 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
2192eec603 Mart*0367       DO ITER=1,IMAX_TICE
16b9a002a2 Mart*0368 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0369        ikey = ITER
                0370 CADJ STORE tsurfLoc = loctape_solve4temp, key = ikey, byte = isbyte
16b9a002a2 Mart*0371 #endif
2192eec603 Mart*0372        DO J=1,sNy
                0373         DO I=1,sNx
6d55e2b45a Mart*0374
112191cc93 Jean*0375 C-    save tsurf from previous iter
                0376          tsurfPrev(I,J) = tsurfLoc(I,J)
2192eec603 Mart*0377          IF ( iceOrNot(I,J) ) THEN
9637aec598 Jean*0378
2192eec603 Mart*0379           t1 = tsurfLoc(I,J)
                0380           t2 = t1*t1
                0381           t3 = t2*t1
                0382           t4 = t2*t2
9637aec598 Jean*0383
6e70381b89 Jean*0384 C--   Calculate the specific humidity in the BL above the snow/ice
112191cc93 Jean*0385           IF ( useMaykutSatVapPoly ) THEN
6e70381b89 Jean*0386 C-    Use the Maykut polynomial
112191cc93 Jean*0387            qhice(I,J)=QS1*(C1*t4+C2*t3 +C3*t2+C4*t1+C5)
                0388            dqh_dTs(I,J) = 0. _d 0
                0389           ELSE
6e70381b89 Jean*0390 C-    Use exponential relation approx., more accurate at low temperatures
2192eec603 Mart*0391 C     log 10 of the sat vap pressure
112191cc93 Jean*0392            mm_log10pi = -aa1 / t1 + aa2
2192eec603 Mart*0393 C     The saturation vapor pressure (SVP) in the surface
                0394 C     boundary layer (BL) above the snow/ice.
112191cc93 Jean*0395 c          mm_pi = TEN **(mm_log10pi)
4fecc0d4b4 Jean*0396 C     The following form does the same, but is faster
112191cc93 Jean*0397            mm_pi = EXP(mm_log10pi*lnTEN)
                0398            qhice(I,J) = bb1*mm_pi/( Ppascals -(1.0 _d 0 - bb1)*mm_pi )
a676451ac2 Jean*0399 C     A constant for SVP derivative w.r.t TICE
112191cc93 Jean*0400 c          cc3t = TEN **(aa1 / t1)
a676451ac2 Jean*0401 C     The following form does the same, but is faster
112191cc93 Jean*0402            cc3t = EXP(aa1 / t1 * lnTEN)
a676451ac2 Jean*0403 C     d(qh)/d(TICE)
112191cc93 Jean*0404            dqh_dTs(I,J) = cc1*cc3t/((cc2-cc3t*Ppascals)**2 *t2)
                0405           ENDIF
9637aec598 Jean*0406
69d8fa39b5 Mart*0407 C     Calculate the flux terms based on the updated tsurfLoc
16b9a002a2 Mart*0408           F_c(I,J)    = effConduct(I,J)*(tempFrz(I,J)-t1)
69d8fa39b5 Mart*0409           F_lh(I,J)   = D1I*UG(I,J)*(qhice(I,J)-AQH(I,J,bi,bj))
840c7fba30 Gael*0410 #ifdef SEAICE_CAP_SUBLIM
a676451ac2 Jean*0411 C     if the latent heat flux implied by tsurfLoc exceeds
f208488703 Mart*0412 C     F_lh_max, cap F_lh and decouple the flux magnitude from tIce (tsurfLoc)
a676451ac2 Jean*0413           IF (F_lh(I,J) .GT. F_lh_max(I,J)) THEN
52ff14d141 Ian *0414              F_lh(I,J)  = F_lh_max(I,J)
6e70381b89 Jean*0415              dqh_dTs(I,J) = ZERO
a676451ac2 Jean*0416           ENDIF
840c7fba30 Gael*0417 #endif /* SEAICE_CAP_SUBLIM */
52ff14d141 Ian *0418
a676451ac2 Jean*0419           F_lwu(I,J) = t4 * D3(I,J)
2192eec603 Mart*0420           F_sens(I,J)= D1 * UG(I,J) * (t1 - atempLoc(I,J))
a676451ac2 Jean*0421           F_ia(I,J) = -lwdownLoc(I,J) -absorbedSW(I,J) + F_lwu(I,J)
6e70381b89 Jean*0422      &              +  F_sens(I,J) + F_lh(I,J)
aa0478ba9b Jean*0423 C     d(F_ia)/d(Tsurf)
                0424           dFia_dTs(I,J) = 4.0 _d 0*D3(I,J)*t3 + D1*UG(I,J)
                0425      &                  + D1I*UG(I,J)*dqh_dTs(I,J)
9637aec598 Jean*0426
                0427 #ifdef SEAICE_DEBUG
a676451ac2 Jean*0428           IF ( (I .EQ. SEAICE_debugPointI) .AND.
52ff14d141 Ian *0429      &         (J .EQ. SEAICE_debugPointJ) ) THEN
2192eec603 Mart*0430            print '(A,i6,4(1x,D24.15))',
                0431      &          'ice-iter qhICE,       ', ITER,qhIce(I,J)
                0432            print '(A,i6,4(1x,D24.15))',
6e70381b89 Jean*0433      &          'ice-iter dFiDTs1 F_ia ', ITER,
                0434      &          dFia_dTs(I,J)+effConduct(I,J), F_ia(I,J)-F_c(I,J)
2192eec603 Mart*0435           ENDIF
300dea8deb Jean*0436 #endif /* SEAICE_DEBUG */
9637aec598 Jean*0437
aa0478ba9b Jean*0438 C-    Update tsurf as solution of : Fc = Fia + d/dT(Fia - Fc) *delta.tsurf
6e70381b89 Jean*0439           tsurfLoc(I,J) = tsurfLoc(I,J)
                0440      &    + ( F_c(I,J)-F_ia(I,J) ) / ( effConduct(I,J)+dFia_dTs(I,J) )
9637aec598 Jean*0441
16b9a002a2 Mart*0442          ENDIF
                0443         ENDDO
                0444        ENDDO
                0445        IF ( useMaykutSatVapPoly ) THEN
aef1621d33 Patr*0446 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0447 CADJ STORE tsurfLoc = loctape_solve4temp, key = ikey, byte = isbyte
16b9a002a2 Mart*0448 #endif
                0449         DO J=1,sNy
                0450          DO I=1,sNx
                0451           tsurfLoc(I,J) = MAX( celsius2K+MIN_TICE, tsurfLoc(I,J) )
                0452 C     If the search leads to tsurfLoc < 50 Kelvin, restart the search at
                0453 C     tsurfLoc = TMELT. Note that one solution to the energy balance problem is
                0454 C     an extremely low temperature - a temperature far below realistic values.
aa0478ba9b Jean*0455 c         IF (tsurfLoc(I,J) .LT. 50.0 _d 0 ) tsurfLoc(I,J) = TMELT
                0456 C   Comments & code above not relevant anymore (from older version, when
                0457 C   trying Maykut-Polynomial & dqh_dTs > 0 ?): commented out
16b9a002a2 Mart*0458          ENDDO
                0459         ENDDO
                0460        ENDIF
                0461 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0462 CADJ STORE tsurfLoc = loctape_solve4temp, key = ikey, byte = isbyte
16b9a002a2 Mart*0463 #endif
                0464        DO J=1,sNy
                0465         DO I=1,sNx
                0466          tsurfLoc(I,J) = MIN( tsurfLoc(I,J), TMELT )
9637aec598 Jean*0467
                0468 #ifdef SEAICE_DEBUG
16b9a002a2 Mart*0469          IF ( (I .EQ. SEAICE_debugPointI) .AND.
                0470      &        (J .EQ. SEAICE_debugPointJ) ) THEN
                0471           print '(A,i6,4(1x,D24.15))',
                0472      &         'ice-iter tsurfLc,|dif|', ITER,
                0473      &         tsurfLoc(I,J),
                0474      &         LOG10(ABS(tsurfLoc(I,J) - tsurfPrev(I,J)))
                0475          ENDIF
300dea8deb Jean*0476 #endif /* SEAICE_DEBUG */
9637aec598 Jean*0477
2192eec603 Mart*0478         ENDDO                   !/* i */
                0479        ENDDO                    !/* j */
                0480       ENDDO                     !/* Iterations */
a676451ac2 Jean*0481 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0482
112191cc93 Jean*0483 #ifdef SEAICE_MODIFY_GROWTH_ADJ
                0484 Cgf no additional dependency through solver, snow, etc.
16b9a002a2 Mart*0485       IF ( SEAICEadjMODE.GE.2 ) THEN
                0486        DO J=1,sNy
                0487         DO I=1,sNx
                0488          IF ( iceOrNot(I,J) ) THEN
f5282c5b03 Gael*0489           CALL ZERO_ADJ_1D( 1, TSURFin(I,J), myThid)
aa0478ba9b Jean*0490           absorbedSW(I,J) = 0.3 _d 0 *SWDOWN(I,J,bi,bj)
112191cc93 Jean*0491           IcePenetSW(I,J)= 0. _d 0
aa0478ba9b Jean*0492          ENDIF
16b9a002a2 Mart*0493         ENDIF
                0494        ENDDO
                0495       ENDDO
                0496       IF ( postSolvTempIter.EQ.2 .OR. SEAICEadjMODE.GE.2 ) THEN
                0497        DO J=1,sNy
                0498         DO I=1,sNx
                0499          IF ( iceOrNot(I,J) ) THEN
f5282c5b03 Gael*0500           t1 = TSURFin(I,J)
112191cc93 Jean*0501 #else /* SEAICE_MODIFY_GROWTH_ADJ */
16b9a002a2 Mart*0502       IF ( postSolvTempIter.EQ.2 ) THEN
                0503        DO J=1,sNy
                0504         DO I=1,sNx
                0505          IF ( iceOrNot(I,J) ) THEN
2192eec603 Mart*0506 C     Recalculate the fluxes based on the (possibly) adjusted TSURF
112191cc93 Jean*0507           t1 = tsurfLoc(I,J)
aa0478ba9b Jean*0508 #endif /* SEAICE_MODIFY_GROWTH_ADJ */
112191cc93 Jean*0509           t2 = t1*t1
                0510           t3 = t2*t1
                0511           t4 = t2*t2
9637aec598 Jean*0512
112191cc93 Jean*0513           IF ( useMaykutSatVapPoly ) THEN
                0514            qhice(I,J)=QS1*(C1*t4+C2*t3 +C3*t2+C4*t1+C5)
                0515           ELSE
2192eec603 Mart*0516 C     log 10 of the sat vap pressure
112191cc93 Jean*0517            mm_log10pi = -aa1 / t1 + aa2
2192eec603 Mart*0518 C     saturation vapor pressure
112191cc93 Jean*0519 c          mm_pi = TEN **(mm_log10pi)
4fecc0d4b4 Jean*0520 C     The following form does the same, but is faster
112191cc93 Jean*0521            mm_pi = EXP(mm_log10pi*lnTEN)
a676451ac2 Jean*0522 C     over ice specific humidity
112191cc93 Jean*0523            qhice(I,J) = bb1*mm_pi/( Ppascals -(1.0 _d 0 - bb1)*mm_pi )
                0524           ENDIF
840c7fba30 Gael*0525           F_c(I,J)  = effConduct(I,J) * (tempFrz(I,J) - t1)
112191cc93 Jean*0526           F_lh(I,J) = D1I * UG(I,J)*(qhice(I,J)-AQH(I,J,bi,bj))
840c7fba30 Gael*0527 #ifdef SEAICE_CAP_SUBLIM
112191cc93 Jean*0528           IF (F_lh(I,J) .GT. F_lh_max(I,J)) THEN
52ff14d141 Ian *0529              F_lh(I,J)  = F_lh_max(I,J)
112191cc93 Jean*0530           ENDIF
840c7fba30 Gael*0531 #endif /* SEAICE_CAP_SUBLIM */
112191cc93 Jean*0532           F_lwu(I,J)  = t4 * D3(I,J)
                0533           F_sens(I,J) = D1 * UG(I,J) * (t1 - atempLoc(I,J))
a676451ac2 Jean*0534 C     The flux between the ice/snow surface and the atmosphere.
112191cc93 Jean*0535           F_ia(I,J) = -lwdownLoc(I,J) -absorbedSW(I,J) + F_lwu(I,J)
                0536      &              +  F_sens(I,J) + F_lh(I,J)
9637aec598 Jean*0537
4dd39c50d9 Mart*0538 C IGF_SIR-b
                0539           IF (-F_c(I,J) .LT. ZERO) THEN
                0540 C note that F_c is flipped sign in this rewrite for some reason
                0541             F_io_net(I,J) = F_c(I,J)
                0542             F_ia_net(I,J) = ZERO
                0543           ELSE
                0544             F_io_net(I,J) = ZERO
                0545             F_ia_net(I,J) = F_ia(I,J)
                0546           ENDIF
                0547 C IGF_SIR-e
                0548
16b9a002a2 Mart*0549          ENDIF
                0550         ENDDO
                0551        ENDDO
                0552       ELSEIF ( postSolvTempIter.EQ.1 ) THEN
                0553        DO J=1,sNy
                0554         DO I=1,sNx
                0555          IF ( iceOrNot(I,J) ) THEN
112191cc93 Jean*0556 C     Update fluxes (consistent with the linearized formulation)
                0557           delTsurf  = tsurfLoc(I,J)-tsurfPrev(I,J)
840c7fba30 Gael*0558           F_c(I,J)  = effConduct(I,J)*(tempFrz(I,J)-tsurfLoc(I,J))
112191cc93 Jean*0559           F_ia(I,J) = F_ia(I,J) + dFia_dTs(I,J)*delTsurf
                0560           F_lh(I,J) = F_lh(I,J)
                0561      &              + D1I*UG(I,J)*dqh_dTs(I,J)*delTsurf
aa0478ba9b Jean*0562
                0563 c        ELSEIF ( postSolvTempIter.EQ.0 ) THEN
112191cc93 Jean*0564 C     Take fluxes from last iteration
aa0478ba9b Jean*0565
a676451ac2 Jean*0566          ENDIF
16b9a002a2 Mart*0567         ENDDO
                0568        ENDDO
                0569       ENDIF
                0570       DO J=1,sNy
                0571        DO I=1,sNx
                0572         IF ( iceOrNot(I,J) ) THEN
a676451ac2 Jean*0573
16b9a002a2 Mart*0574 C     Save updated tsurf and finalize the flux terms
                0575          TSURFout(I,J) = tsurfLoc(I,J)
a676451ac2 Jean*0576 C     Fresh water flux (kg/m^2/s) from latent heat of sublimation.
                0577 C     F_lh is positive upward (sea ice looses heat) and FWsublim
                0578 C     is also positive upward (atmosphere gains freshwater)
                0579          FWsublim(I,J) = F_lh(I,J)/lhSublim
3a3bf6419a Gael*0580
a676451ac2 Jean*0581 #ifdef SEAICE_DEBUG
                0582          IF ( (I .EQ. SEAICE_debugPointI) .AND.
52ff14d141 Ian *0583      &        (J .EQ. SEAICE_debugPointJ) ) THEN
2192eec603 Mart*0584           print '(A)','----------------------------------------'
                0585           print '(A,i6)','ibi complete ', myIter
                0586           print '(A,4(1x,D24.15))',
                0587      &         'ibi T(SURF, surfLoc,atmos) ',
f5282c5b03 Gael*0588      &         TSURFout(I,J), tsurfLoc(I,J),atempLoc(I,J)
2192eec603 Mart*0589           print '(A,4(1x,D24.15))',
                0590      &         'ibi LWL                    ', lwdownLoc(I,J)
                0591           print '(A,4(1x,D24.15))',
                0592      &         'ibi QSW(Total, Penetrating)',
a676451ac2 Jean*0593      &         SWDOWN(I,J,bi,bj), IcePenetSW(I,J)
2192eec603 Mart*0594           print '(A,4(1x,D24.15))',
                0595      &         'ibi qh(ATM ICE)            ',
                0596      &         AQH(I,J,bi,bj),qhice(I,J)
52ff14d141 Ian *0597          print '(A,4(1x,D24.15))',
                0598      &         'ibi F(lwd,swi,lwu)         ',
a676451ac2 Jean*0599      &         -lwdownLoc(I,J), -absorbedSW(I,J), F_lwu(I,J)
52ff14d141 Ian *0600          print '(A,4(1x,D24.15))',
                0601      &         'ibi F(c,lh,sens)           ',
                0602      &         F_c(I,J), F_lh(I,J), F_sens(I,J)
840c7fba30 Gael*0603 #ifdef SEAICE_CAP_SUBLIM
52ff14d141 Ian *0604          IF (F_lh_max(I,J) .GT. ZERO) THEN
                0605              print '(A,4(1x,D24.15))',
                0606      &         'ibi F_lh_max,  F_lh/lhmax) ',
                0607      &         F_lh_max(I,J), F_lh(I,J)/ F_lh_max(I,J)
1043a55cc1 Jean*0608          ELSE
52ff14d141 Ian *0609              print '(A,4(1x,D24.15))',
                0610      &         'ibi F_lh_max = ZERO! '
                0611          ENDIF
                0612          print '(A,4(1x,D24.15))',
                0613      &         'ibi FWsub, FWsubm*dT/rhoI  ',
                0614      &          FWsublim(I,J),
                0615      &          FWsublim(I,J)*SEAICE_deltaTtherm/SEAICE_rhoICE
840c7fba30 Gael*0616 #endif /* SEAICE_CAP_SUBLIM */
2192eec603 Mart*0617           print '(A,4(1x,D24.15))',
                0618      &         'ibi F_ia, F_ia_net, F_c    ',
4dd39c50d9 Mart*0619      &         F_ia(I,J), F_ia_net(I,J), F_c(I,J)
2192eec603 Mart*0620           print '(A)','----------------------------------------'
                0621          ENDIF
300dea8deb Jean*0622 #endif /* SEAICE_DEBUG */
4fecc0d4b4 Jean*0623
2192eec603 Mart*0624         ENDIF                   !/* iceOrNot */
                0625        ENDDO                    !/* i */
9637aec598 Jean*0626       ENDDO                     !/* j */
                0627
1043a55cc1 Jean*0628 #endif /* ALLOW_ATM_TEMP && ALLOW_DOWNWARD_RADIATION */
                0629       RETURN
9637aec598 Jean*0630       END