/[MITgcm]/MITgcm/pkg/thsice/thsice_solve4temp.F
ViewVC logotype

Diff of /MITgcm/pkg/thsice/thsice_solve4temp.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.15 by heimbach, Tue Apr 17 23:42:33 2007 UTC revision 1.24 by heimbach, Sat Oct 16 19:22:34 2010 UTC
# Line 7  CBOP Line 7  CBOP
7  C     !ROUTINE: THSICE_SOLVE4TEMP  C     !ROUTINE: THSICE_SOLVE4TEMP
8  C     !INTERFACE:  C     !INTERFACE:
9        SUBROUTINE THSICE_SOLVE4TEMP(        SUBROUTINE THSICE_SOLVE4TEMP(
10       I                  bi, bj, siLo, siHi, sjLo, sjHi,       I                  bi, bj,
11       I                  iMin,iMax, jMin,jMax, dBugFlag,       I                  iMin,iMax, jMin,jMax, dBugFlag,
12       I                  useBulkForce, useEXF,       I                  useBulkForce, useEXF,
13       I                  iceMask, hIce, hSnow, tFrz, flxExSW,       I                  iceMask, hIce, hSnow, tFrz, flxExSW,
14       U                  flxSW, tSrf, qIc1, qIc2,       U                  flxSW, tSrf, qIc1, qIc2,
# Line 27  C ADAPTED FROM: Line 27  C ADAPTED FROM:
27  C LANL CICE.v2.0.2  C LANL CICE.v2.0.2
28  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
29  C.. thermodynamics (vertical physics) based on M. Winton 3-layer model  C.. thermodynamics (vertical physics) based on M. Winton 3-layer model
30  C.. See Bitz, C. M. and W. H. Lipscomb, 1999:  "An energy-conserving  C.. See Bitz, C. M. and W. H. Lipscomb, 1999:  An energy-conserving
31  C..       thermodynamic sea ice model for climate study."  J. Geophys.  C..       thermodynamic sea ice model for climate study.
32  C..       Res., 104, 15669 - 15677.  C..       J. Geophys. Res., 104, 15669 - 15677.
33  C..     Winton, M., 1999:  "A reformulated three-layer sea ice model."  C..     Winton, M., 1999:  "A reformulated three-layer sea ice model."
34  C..       Submitted to J. Atmos. Ocean. Technol.  C..       Submitted to J. Atmos. Ocean. Technol.
35  C.. authors Elizabeth C. Hunke and William Lipscomb  C.. authors Elizabeth C. Hunke and William Lipscomb
# Line 44  C     !USES: Line 44  C     !USES:
44    
45  C     == Global variables ===  C     == Global variables ===
46  #include "EEPARAMS.h"  #include "EEPARAMS.h"
47    #include "SIZE.h"
48  #include "THSICE_SIZE.h"  #include "THSICE_SIZE.h"
49  #include "THSICE_PARAMS.h"  #include "THSICE_PARAMS.h"
50  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
 # include "SIZE.h"  
51  # include "tamc.h"  # include "tamc.h"
52  # include "tamc_keys.h"  # include "tamc_keys.h"
53  #endif  #endif
54    
55  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
56  C     == Routine Arguments ==  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  
57  C     bi,bj       :: tile indices  C     bi,bj       :: tile indices
58  C     iMin,iMax   :: computation domain: 1rst index range  C     iMin,iMax   :: computation domain: 1rst index range
59  C     jMin,jMax   :: computation domain: 2nd  index range  C     jMin,jMax   :: computation domain: 2nd  index range
# Line 63  C     dBugFlag    :: allow to print debu Line 61  C     dBugFlag    :: allow to print debu
61  C     useBulkForce:: use surf. fluxes from bulk-forcing external S/R  C     useBulkForce:: use surf. fluxes from bulk-forcing external S/R
62  C     useEXF      :: use surf. fluxes from exf          external S/R  C     useEXF      :: use surf. fluxes from exf          external S/R
63  C---  Input:  C---  Input:
64  C         iceMask :: sea-ice fractional mask [0-1]  C     iceMask     :: sea-ice fractional mask [0-1]
65  C  hIce    (hi)   :: ice height [m]  C     hIce        :: ice height [m]
66  C  hSnow   (hs)   :: snow height [m]  C     hSnow       :: snow height [m]
67  C  tFrz    (Tf)   :: sea-water freezing temperature [oC] (function of S)  C     tFrz        :: sea-water freezing temperature [oC] (function of S)
68  C  flxExSW  (=)   :: surf. heat flux (+=down) except SW, function of surf. temp Ts:  C     flxExSW     :: surf. heat flux (+=down) except SW, function of surf. temp Ts:
69  C                    0: Flx(Ts=0) ; 1: Flx(Ts=Ts^n) ; 2: d.Flx/dTs(Ts=Ts^n)  C                    0: Flx(Ts=0) ; 1: Flx(Ts=Ts^n) ; 2: d.Flx/dTs(Ts=Ts^n)
70  C---  Modified (input&output):  C---  Modified (input&output):
71  C  flxSW  (netSW) :: net Short-Wave flux (+=down) [W/m2]: input= at surface  C     flxSW       :: net Short-Wave flux (+=down) [W/m2]: input= at surface
72  C           (=)   ::                 output= below sea-ice, into the ocean  C                 ::                 output= below sea-ice, into the ocean
73  C  tSrf    (Tsf)  :: surface (ice or snow) temperature  C     tSrf        :: surface (ice or snow) temperature
74  C  qIc1   (qicen) :: ice enthalpy (J/kg), 1rst level  C     qIc1        :: ice enthalpy (J/kg), 1rst level
75  C  qIc2   (qicen) :: ice enthalpy (J/kg), 2nd level  C     qIc2        :: ice enthalpy (J/kg), 2nd level
76  C---  Output  C---  Output
77  C  tIc1    (Tice) :: temperature of ice layer 1 [oC]  C     tIc1        :: temperature of ice layer 1 [oC]
78  C  tIc2    (Tice) :: temperature of ice layer 2 [oC]  C     tIc2        :: temperature of ice layer 2 [oC]
79  C  dTsrf   (dTsf) :: surf. temp adjusment: Ts^n+1 - Ts^n  C     dTsrf       :: surf. temp adjusment: Ts^n+1 - Ts^n
80  C  sHeat(sHeating):: surf heating flux left to melt snow or ice (= Atmos-conduction)  C     sHeat       :: surf heating flux left to melt snow or ice (= Atmos-conduction)
81  C  flxCnB   (=)   :: heat flux conducted through the ice to bottom surface  C     flxCnB      :: heat flux conducted through the ice to bottom surface
82  C  flxAtm   (=)   :: net flux of energy from the atmosphere [W/m2] (+=down)  C     flxAtm      :: net flux of energy from the atmosphere [W/m2] (+=down)
83  C                    without snow precip. (energy=0 for liquid water at 0.oC)  C                    without snow precip. (energy=0 for liquid water at 0.oC)
84  C  evpAtm   (=)   :: evaporation to the atmosphere [kg/m2/s] (>0 if evaporate)  C     evpAtm      :: evaporation to the atmosphere [kg/m2/s] (>0 if evaporate)
85  C---  Input:  C---  Input:
86  C     myTime      :: current Time of simulation [s]  C     myTime      :: current Time of simulation [s]
87  C     myIter      :: current Iteration number in simulation  C     myIter      :: current Iteration number in simulation
88  C     myThid      :: my Thread Id number  C     myThid      :: my Thread Id number
       INTEGER siLo, siHi, sjLo, sjHi  
89        INTEGER bi,bj        INTEGER bi,bj
90        INTEGER iMin, iMax        INTEGER iMin, iMax
91        INTEGER jMin, jMax        INTEGER jMin, jMax
92        LOGICAL dBugFlag        LOGICAL dBugFlag
93        LOGICAL useBulkForce        LOGICAL useBulkForce
94        LOGICAL useEXF        LOGICAL useEXF
95        _RL iceMask(siLo:siHi,sjLo:sjHi)        _RL iceMask(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
96        _RL hIce   (siLo:siHi,sjLo:sjHi)        _RL hIce   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
97        _RL hSnow  (siLo:siHi,sjLo:sjHi)        _RL hSnow  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98        _RL tFrz   (siLo:siHi,sjLo:sjHi)        _RL tFrz   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
99        _RL flxExSW(iMin:iMax,jMin:jMax,0:2)        _RL flxExSW(iMin:iMax,jMin:jMax,0:2)
100        _RL flxSW  (siLo:siHi,sjLo:sjHi)        _RL flxSW  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
101        _RL tSrf   (siLo:siHi,sjLo:sjHi)        _RL tSrf   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102        _RL qIc1   (siLo:siHi,sjLo:sjHi)        _RL qIc1   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103        _RL qIc2   (siLo:siHi,sjLo:sjHi)        _RL qIc2   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
104        _RL tIc1   (siLo:siHi,sjLo:sjHi)        _RL tIc1   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105        _RL tIc2   (siLo:siHi,sjLo:sjHi)        _RL tIc2   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
 c     _RL dTsrf  (siLo:siHi,sjLo:sjHi)  
106        _RL dTsrf  (iMin:iMax,jMin:jMax)        _RL dTsrf  (iMin:iMax,jMin:jMax)
107        _RL sHeat  (siLo:siHi,sjLo:sjHi)        _RL sHeat  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
108        _RL flxCnB (siLo:siHi,sjLo:sjHi)        _RL flxCnB (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
109        _RL flxAtm (siLo:siHi,sjLo:sjHi)        _RL flxAtm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
110        _RL evpAtm (siLo:siHi,sjLo:sjHi)        _RL evpAtm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
111        _RL  myTime        _RL  myTime
112        INTEGER myIter        INTEGER myIter
113        INTEGER myThid        INTEGER myThid
# Line 119  CEOP Line 115  CEOP
115    
116  #ifdef ALLOW_THSICE  #ifdef ALLOW_THSICE
117  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
 C---  local copy of input/output argument list variables (see description above)  
 c     _RL flxExcSw(0:2)  
       _RL Tf  
       _RL hi  
       _RL hs  
       _RL netSW  
       _RL Tsf  
       _RL qicen(nlyr)  
       _RL Tice (nlyr)  
 c     _RL sHeating  
 c     _RL flxCnB  
       _RL dTsf  
 c     _RL flxAtm  
 c     _RL evpAtm  
   
118  C     == Local Variables ==  C     == Local Variables ==
119        INTEGER i,j  C     useBlkFlx    :: use some bulk-formulae to compute fluxes over sea-ice
120        INTEGER k, iterMax  C                  :: (otherwise, fluxes are passed as argument from AIM)
121    C     iterate4Tsf  :: to stop to iterate when all icy grid pts Tsf did converged
122        _RL  frsnow        ! fractional snow cover  C     iceFlag      :: True= do iterate for Surf.Temp ; False= do nothing
123        _RL  fswpen        ! SW penetrating beneath surface (W m-2)  C     frsnow       :: fractional snow cover
124        _RL  fswdn         ! SW absorbed at surface (W m-2)  C     fswpen       :: SW penetrating beneath surface (W m-2)
125        _RL  fswint        ! SW absorbed in ice (W m-2)  C     fswdn        :: SW absorbed at surface (W m-2)
126        _RL  fswocn        ! SW passed through ice to ocean (W m-2)  C     fswint       :: SW absorbed in ice (W m-2)
127        _RL  flxExceptSw   ! net surface heat flux, except short-wave (W/m2)  C     fswocn       :: SW passed through ice to ocean (W m-2)
128  C     evap           :: evaporation over snow/ice [kg/m2/s] (>0 if evaporate)  C     Tsf          :: surface (ice or snow) temperature   (local copy of tSrf)
129  C     dEvdT          :: derivative of evap. with respect to Tsf [kg/m2/s/K]  C     flx0exSW     :: net surface heat flux over melting snow/ice, except short-wave.
130        _RL  evap, dEvdT  C     flxTexSW     :: net surface heat flux, except short-wave (W/m2)
131        _RL  flx0          ! net surf heat flux, from Atmos. to sea-ice (W m-2)  C     dFlxdT       :: deriv of flxNet wrt Tsf (W m-2 deg-1)
132        _RL  fct           ! heat conducted to top surface  C     evap0        :: evaporation over melting snow/ice [kg/m2/s] (>0 if evaporate)
133        _RL  df0dT         ! deriv of flx0 wrt Tsf (W m-2 deg-1)  C     evapT        :: evaporation over snow/ice [kg/m2/s] (>0 if evaporate)
134        _RL  k12, k32      ! thermal conductivity terms  C     dEvdT        :: derivative of evap. with respect to Tsf [kg/m2/s/K]
135        _RL  a10, b10      ! coefficients in quadratic eqn for T1  C     flxNet       :: net surf heat flux (+=down), from Atmos. to sea-ice (W m-2)
136        _RL  a1, b1, c1    ! coefficients in quadratic eqn for T1  C     netSW        :: net Short-Wave flux at surface (+=down) [W/m2]
137  c     _RL  Tsf_start     ! old value of Tsf  C     fct          :: heat conducted to top surface
138        _RL  dt            ! timestep  C     k12, k32     :: thermal conductivity terms
139        INTEGER iceornot  C     a10,b10,c10  :: coefficients in quadratic eqn for T1
140    C     a1, b1, c1   :: coefficients in quadratic eqn for T1
141    C     dt           :: timestep
142        LOGICAL useBlkFlx        LOGICAL useBlkFlx
143          LOGICAL iterate4Tsf
144          LOGICAL iceFlag(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
145          INTEGER i,j
146          INTEGER k, iterMax, ii, jj, icount
147          _RL  frsnow
148          _RL  fswpen
149          _RL  fswdn
150          _RL  fswint
151          _RL  fswocn
152          _RL  Tsf     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
153          _RL  flx0exSW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
154          _RL  flxTexSW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
155          _RL  dFlxdT  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
156          _RL  evap0   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
157          _RL  evapT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
158          _RL  dEvdT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
159          _RL  flxNet
160          _RL  fct
161          _RL  k12     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
162          _RL  k32
163          _RL  a10     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
164          _RL  b10     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
165          _RL  c10     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
166          _RL  a1, b1, c1
167          _RL  dt
168          _RL  recip_dhSnowLin
169    #ifdef ALLOW_DBUG_THSICE
170          _RL  netSW
171    #endif
172    
173  C-    define grid-point location where to print debugging values  C-    Define grid-point location where to print debugging values
174  #include "THSICE_DEBUG.h"  #include "THSICE_DEBUG.h"
175    
176   1010 FORMAT(A,I3,3F11.6)   1010 FORMAT(A,I3,3F11.6)
177   1020 FORMAT(A,1P4E14.6)   1020 FORMAT(A,1P4E14.6)
178    
179  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|
180    
181  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
182        act1 = bi - myBxLo(myThid)  c     act1 = bi - myBxLo(myThid)
183        max1 = myBxHi(myThid) - myBxLo(myThid) + 1  c     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
184        act2 = bj - myByLo(myThid)  c     act2 = bj - myByLo(myThid)
185        max2 = myByHi(myThid) - myByLo(myThid) + 1  c     max2 = myByHi(myThid) - myByLo(myThid) + 1
186        act3 = myThid - 1  c     act3 = myThid - 1
187        max3 = nTx*nTy  c     max3 = nTx*nTy
188        act4 = ikey_dynamics - 1  c     act4 = ikey_dynamics - 1
189          iicekey = (act1 + 1) + act2*max1
190         &                     + act3*max1*max2
191         &                     + act4*max1*max2*max3
192  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
193    
194        useBlkFlx = useEXF .OR. useBulkForce  #ifdef ALLOW_AUTODIFF_TAMC
195    CADJ STORE flxsw(:,:) = comlev1_bibj,key=iicekey,byte=isbyte
196          DO j = jMin, jMax
197           DO i = iMin, iMax
198            tic1(i,j) = 0. _d 0
199            tic2(i,j) = 0. _d 0
200           END DO
201          END DO
202    #endif
203    
204          useBlkFlx = useEXF .OR. useBulkForce
205        dt  = thSIce_dtTemp        dt  = thSIce_dtTemp
206          IF ( dhSnowLin.GT.0. _d 0 ) THEN
207            recip_dhSnowLin = 1. _d 0 / dhSnowLin
208          ELSE
209            recip_dhSnowLin = 0. _d 0
210          ENDIF
211    
212          iterate4Tsf = .FALSE.
213          icount = 0
214    C    
215        DO j = jMin, jMax        DO j = jMin, jMax
216         DO i = iMin, iMax         DO i = iMin, iMax
217  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
218            ikey_1 = i  c         ikey_1 = i
219       &         + sNx*(j-1)  c    &         + sNx*(j-1)
220       &         + sNx*sNy*act1  c    &         + sNx*sNy*act1
221       &         + sNx*sNy*max1*act2  c    &         + sNx*sNy*max1*act2
222       &         + sNx*sNy*max1*max2*act3  c    &         + sNx*sNy*max1*max2*act3
223       &         + sNx*sNy*max1*max2*max3*act4  c    &         + sNx*sNy*max1*max2*max3*act4
224  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
 C--  
225  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
226  CADJ STORE  devdt        = comlev1_thsice_1, key=ikey_1  cCADJ STORE  devdt        = comlev1_thsice_1, key=ikey_1
227  CADJ STORE  df0dt        = comlev1_thsice_1, key=ikey_1  cCADJ STORE  dFlxdT        = comlev1_thsice_1, key=ikey_1
228  CADJ STORE  flxexceptsw  = comlev1_thsice_1, key=ikey_1  cCADJ STORE  flxexceptsw  = comlev1_thsice_1, key=ikey_1
229  CADJ STORE  flxsw(i,j)   = comlev1_thsice_1, key=ikey_1  cCADJ STORE  flxsw(i,j)   = comlev1_thsice_1, key=ikey_1
230  CADJ STORE  qic1(i,j)    = comlev1_thsice_1, key=ikey_1  cCADJ STORE  qic1(i,j)    = comlev1_thsice_1, key=ikey_1
231  CADJ STORE  qic2(i,j)    = comlev1_thsice_1, key=ikey_1  cCADJ STORE  qic2(i,j)    = comlev1_thsice_1, key=ikey_1
232  CADJ STORE  tsrf(i,j)    = comlev1_thsice_1, key=ikey_1  cCADJ STORE  tsrf(i,j)    = comlev1_thsice_1, key=ikey_1
233  #endif  #endif
234          IF ( iceMask(i,j).GT.0. _d 0) THEN          IF ( iceMask(i,j).GT.0. _d 0) THEN
235            hi      = hIce(i,j)            iterate4Tsf  = .TRUE.
236            hs      = hSnow(i,j)            iceFlag(i,j) = .TRUE.
           Tf      = tFrz(i,j)  
           netSW   = flxSW(i,j)  
           Tsf     = tSrf(i,j)  
           qicen(1)= qIc1(i,j)  
           qicen(2)= qIc2(i,j)  
 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  
237  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
238            IF ( dBug(i,j,bi,bj) ) WRITE(6,'(A,2I4,2I2)')            IF ( dBug(i,j,bi,bj) ) WRITE(6,'(A,2I4,2I2)')
239       &         'ThSI_SOLVE4T: i,j=',i,j,bi,bj       &         'ThSI_SOLVE4T: i,j=',i,j,bi,bj
240  #endif  #endif
241            IF ( hi.LT.himin ) THEN            IF ( hIce(i,j).LT.hIceMin ) THEN
242  C If hi < himin, melt the ice.  C     if hi < hIceMin, melt the ice.
243               STOP 'THSICE_SOLVE4TEMP: should not enter if hi<himin'  C     keep the position of this problem but do the stop later
244               ii = i
245               jj = j
246               icount = icount + 1
247              ENDIF
248    
249    C--   Fractional snow cover:
250    C     assume a linear distribution of snow thickness, with dhSnowLin slope,
251    C      from hs-dhSnowLin to hs+dhSnowLin if full ice & snow cover.
252    C     frsnow = fraction of snow over the ice-covered part of the grid cell
253              IF ( hSnow(i,j) .GT. iceMask(i,j)*dhSnowLin ) THEN
254               frsnow = 1. _d 0
255              ELSE
256               frsnow = hSnow(i,j)*recip_dhSnowLin/iceMask(i,j)
257               IF ( frsnow.GT.0. _d 0 ) frsnow = SQRT(frsnow)
258            ENDIF            ENDIF
259    
260  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C--   Compute SW flux absorbed at surface and penetrating to layer 1.
261              fswpen = flxSW(i,j) * (1. _d 0 - frsnow) * i0swFrac
262              fswocn = fswpen * exp(-ksolar*hIce(i,j))
263              fswint = fswpen - fswocn
264              fswdn  = flxSW(i,j) - fswpen
265    
266    C     Initialise Atmospheric surf. heat flux with net SW flux at the surface
267              flxAtm(i,j) = flxSW(i,j)
268    C     SW flux at sea-ice base left to the ocean
269              flxSW(i,j) = fswocn
270    C     Initialise surface Heating with SW contribution
271              sHeat(i,j) = fswdn
272    
273    C--   Compute conductivity terms at layer interfaces.
274              k12(i,j) = 4. _d 0*kIce*kSnow
275         &             / (kSnow*hIce(i,j) + 4. _d 0*kIce*hSnow(i,j))
276              k32      = 2. _d 0*kIce  / hIce(i,j)
277    
278    C--   Compute ice temperatures
279              a1 = cpIce
280              b1 = qIc1(i,j) + (cpWater-cpIce )*Tmlt1 - Lfresh
281              c1 = Lfresh * Tmlt1
282              tIc1(i,j) = 0.5 _d 0 *(-b1 - SQRT(b1*b1-4. _d 0*a1*c1))/a1
283              tIc2(i,j) = (Lfresh-qIc2(i,j)) / cpIce
284    
 C fractional snow cover  
       frsnow = 0. _d 0  
       IF (hs .GT. 0. _d 0) frsnow = 1. _d 0  
   
 C Compute SW flux absorbed at surface and penetrating to layer 1.  
       fswpen  = netSW * (1. _d 0 - frsnow) * i0  
       fswocn = fswpen * exp(-ksolar*hi)  
       fswint = fswpen - fswocn  
   
       fswdn = netSW - fswpen  
   
 C Compute conductivity terms at layer interfaces.  
   
       k12 = 4. _d 0*kice*ksnow / (ksnow*hi + 4. _d 0*kice*hs)  
       k32 = 2. _d 0*kice  / hi  
   
 C compute ice temperatures  
       a1 = cpice  
       b1 = qicen(1) + (cpwater-cpice )*Tmlt1 - Lfresh  
       c1 = Lfresh * Tmlt1  
       Tice(1) = 0.5 _d 0 *(-b1 - SQRT(b1*b1-4. _d 0*a1*c1))/a1  
       Tice(2) = (Lfresh-qicen(2)) / cpice  
   
       IF (Tice(1).GT.0. _d 0 ) THEN  
        WRITE (standardMessageUnit,'(A,I12,1PE14.6)')  
      &       ' BBerr: Tice(1) > 0 ; it=', myIter, qicen(1)  
        WRITE (standardMessageUnit,'(A,4I5,2F11.4)')  
      &      ' BBerr: i,j,bi,bj,Tice = ',i,j,bi,bj,Tice  
       ENDIF  
       IF ( Tice(2).GT.0. _d 0) THEN  
        WRITE (standardMessageUnit,'(A,I12,1PE14.6)')  
      &       ' BBerr: Tice(2) > 0 ; it=', myIter, qicen(2)  
        WRITE (standardMessageUnit,'(A,4I5,2F11.4)')  
      &      ' BBerr: i,j,bi,bj,Tice = ',i,j,bi,bj,Tice  
       ENDIF  
285  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
286        IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)            IF (tIc1(i,j).GT.0. _d 0 ) THEN
287       &  'ThSI_SOLVE4T: k, Ts, Tice=',0,Tsf,Tice             WRITE (standardMessageUnit,'(A,I12,1PE14.6)')
288         &       ' BBerr: Tice(1) > 0 ; it=', myIter, qIc1(i,j)
289               WRITE (standardMessageUnit,'(A,4I5,2F11.4)')
290         &      ' BBerr: i,j,bi,bj,Tice = ',i,j,bi,bj,tIc1(i,j),tIc2(i,j)
291              ENDIF
292              IF ( tIc2(i,j).GT.0. _d 0) THEN
293               WRITE (standardMessageUnit,'(A,I12,1PE14.6)')
294         &       ' BBerr: Tice(2) > 0 ; it=', myIter, qIc2(i,j)
295               WRITE (standardMessageUnit,'(A,4I5,2F11.4)')
296         &      ' BBerr: i,j,bi,bj,Tice = ',i,j,bi,bj,tIc1(i,j),tIc2(i,j)
297              ENDIF
298              IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)
299         &      'ThSI_SOLVE4T: k, Ts, Tice=',0,tSrf(i,j),tIc1(i,j),tIc2(i,j)
300  #endif  #endif
301    
302  C Compute coefficients used in quadratic formula.  C--   Compute coefficients used in quadratic formula.
303    
304        a10 = rhoi*cpice *hi/(2. _d 0*dt) +            a10(i,j) = rhoi*cpIce *hIce(i,j)/(2. _d 0*dt) +
305       &      k32 * (4. _d 0*dt*k32 + rhoi*cpice *hi)       &          k32*( 4. _d 0*dt*k32 + rhoi*cpIce *hIce(i,j) )
306       &       / (6. _d 0*dt*k32 + rhoi*cpice *hi)       &           / ( 6. _d 0*dt*k32 + rhoi*cpIce *hIce(i,j) )
307        b10 = -hi*            b10(i,j) = -hIce(i,j)*
308       &      (rhoi*cpice*Tice(1)+rhoi*Lfresh*Tmlt1/Tice(1))       &          ( rhoi*cpIce*tIc1(i,j) + rhoi*Lfresh*Tmlt1/tIc1(i,j) )
309       &       /(2. _d 0*dt)       &           /(2. _d 0*dt)
310       &      - k32 * (4. _d 0*dt*k32*Tf+rhoi*cpice *hi*Tice(2))       &        - k32*( 4. _d 0*dt*k32*tFrz(i,j)
311       &       / (6. _d 0*dt*k32 + rhoi*cpice *hi) - fswint       &               +rhoi*cpIce*hIce(i,j)*tIc2(i,j) )
312        c1 = rhoi*Lfresh*hi*Tmlt1 / (2. _d 0*dt)       &           / ( 6. _d 0*dt*k32 + rhoi*cpIce *hIce(i,j) )
313         &        - fswint
314  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|            c10(i,j) = rhoi*Lfresh*hIce(i,j)*Tmlt1 / (2. _d 0*dt)
 C Compute new surface and internal temperatures; iterate until  
 C Tsfc converges.  
315    
316        IF ( useBlkFlx ) THEN          ELSE
317          iterMax = nitMaxTsf            iceFlag(i,j) = .FALSE.
318        ELSE          ENDIF
319          iterMax = 1         ENDDO
320          ENDDO
321          IF ( icount .gt. 0 ) THEN
322           WRITE (standardMessageUnit,'(A,I5,A)')
323         &      'THSICE_SOLVE4TEMP: there are ',icount,
324         &      ' case(s) where hIce<hIceMin;'
325           WRITE (standardMessageUnit,'(A,I3,A1,I3,A)')
326         &      'THSICE_SOLVE4TEMP: the last one was at (',ii,',',jj,
327         &      ') with hIce = ', hIce(ii,jj)
328           STOP 'THSICE_SOLVE4TEMP: should not enter if hIce<hIceMin'
329        ENDIF        ENDIF
       dTsf = Terrmax  
330    
331  C ----- begin iteration  -----  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|
       DO k = 1,iterMax  
332    
333  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
334         ikey_3 = (ikey_1-1)*MaxTsf + k  CADJ STORE devdt    = comlev1_bibj,key=iicekey,byte=isbyte
335    CADJ STORE tsf      = comlev1_bibj,key=iicekey,byte=isbyte
336  #endif  #endif
337    
338    C--   Get surface fluxes over melting surface
339  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
340  CADJ STORE tsf          = comlev1_thsice_3, key=ikey_3        IF ( useBlkFlx ) THEN
341  CADJ STORE dtsf         = comlev1_thsice_3, key=ikey_3  #else
342  CADJ STORE df0dt        = comlev1_thsice_3, key=ikey_3        IF ( useBlkFlx .AND. iterate4Tsf  ) THEN
 CADJ STORE flxexceptsw  = comlev1_thsice_3, key=ikey_3  
343  #endif  #endif
344         IF ( ABS(dTsf).GE.Terrmax ) THEN          DO j = jMin, jMax
345             DO i = iMin, iMax
346               Tsf(i,j) = 0.
347             ENDDO
348            ENDDO
349            IF ( useEXF ) THEN
350               CALL THSICE_GET_EXF(   bi, bj,
351         I                            iMin, iMax, jMin, jMax,
352         I                            iceFlag, hSnow, Tsf,
353         O                            flx0exSW, dFlxdT, evap0, dEvdT,
354         I                            myTime, myIter, myThid )
355    C-    could add this "ifdef" to hide THSICE_GET_BULKF from TAF
356    #ifdef ALLOW_BULK_FORCE
357            ELSEIF ( useBulkForce ) THEN
358               CALL THSICE_GET_BULKF( bi, bj,
359         I                            iMin, iMax, jMin, jMax,
360         I                            iceFlag, hSnow, Tsf,
361         O                            flx0exSW, dFlxdT, evap0, dEvdT,
362         I                            myTime, myIter, myThid )
363    #endif /* ALLOW_BULK_FORCE */
364            ENDIF
365          ENDIF
366    
367  C Save temperatures at start of iteration.  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|
368  c        Tsf_start = Tsf  C---  Compute new surface and internal temperatures; iterate until
369    C     Tsfc converges.
370          DO j = jMin, jMax
371            DO i = iMin, iMax
372              Tsf(i,j)  = tSrf(i,j)
373              dTsrf(i,j) = Terrmax
374            ENDDO
375          ENDDO
376          IF ( useBlkFlx ) THEN
377            iterMax = nitMaxTsf
378          ELSE
379            iterMax = 1
380          ENDIF
381    
         IF ( useBlkFlx ) THEN  
382  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
383  CADJ STORE tsf          = comlev1_thsice_3, key=ikey_3  CADJ STORE devdt    = comlev1_bibj,key=iicekey,byte=isbyte
384  #endif  CADJ STORE dflxdt   = comlev1_bibj,key=iicekey,byte=isbyte
385  C Compute top surface flux.  CADJ STORE flx0exsw = comlev1_bibj,key=iicekey,byte=isbyte
386           IF (hs.GT.3. _d -1) THEN  CADJ STORE flxtexsw = comlev1_bibj,key=iicekey,byte=isbyte
               iceornot=2  
          ELSE  
               iceornot=1  
          ENDIF  
          IF ( useBulkForce ) THEN  
           CALL THSICE_GET_BULKF(  
      I                         iceornot, Tsf,  
      O                         flxExceptSw, df0dT, evap, dEvdT,  
      I                         i,j,bi,bj,myThid )  
          ELSEIF ( useEXF ) THEN  
           CALL THSICE_GET_EXF  (  
      I                         iceornot, Tsf,  
      O                         flxExceptSw, df0dT, evap, dEvdT,  
      I                         i,j,bi,bj,myThid )  
          ENDIF  
         ELSE  
          flxExceptSw = flxExSW(i,j,1)  
          df0dT = flxExSW(i,j,2)  
         ENDIF  
         flx0 = fswdn + flxExceptSw  
 #ifdef ALLOW_DBUG_THSICE  
       IF ( dBug(i,j,bi,bj) ) WRITE(6,1020)  
      &  'ThSI_SOLVE4T: flx0,df0dT,k12,D=', flx0,df0dT,k12,k12-df0dT  
387  #endif  #endif
388    
389  C Compute new top layer and surface temperatures.  C ----- begin iteration  -----
390  C If Tsfc is computed to be > 0 C, fix Tsfc = 0 and recompute T1        DO k = 1,iterMax
391  C with different coefficients.  #ifndef ALLOW_AUTODIFF_TAMC
392           IF ( iterate4Tsf ) THEN
393           iterate4Tsf = .FALSE.
394    #endif
395    
396           IF ( useBlkFlx ) THEN
397    C--   Compute top surface flux.
398             IF ( useEXF ) THEN
399               CALL THSICE_GET_EXF(   bi, bj,
400         I                            iMin, iMax, jMin, jMax,
401         I                            iceFlag, hSnow, Tsf,
402         O                            flxTexSW, dFlxdT, evapT, dEvdT,
403         I                            myTime, myIter, myThid )
404    C-    could add this "ifdef" to hide THSICE_GET_BULKF from TAF
405    #ifdef ALLOW_BULK_FORCE
406             ELSEIF ( useBulkForce ) THEN
407               CALL THSICE_GET_BULKF( bi, bj,
408         I                            iMin, iMax, jMin, jMax,
409         I                            iceFlag, hSnow, Tsf,
410         O                            flxTexSW, dFlxdT, evapT, dEvdT,
411         I                            myTime, myIter, myThid )
412    #endif /* ALLOW_BULK_FORCE */
413             ENDIF
414           ELSE
415             DO j = jMin, jMax
416              DO i = iMin, iMax
417               IF ( iceFlag(i,j) ) THEN
418                 flxTexSW(i,j) = flxExSW(i,j,1)
419                 dFlxdT(i,j) = flxExSW(i,j,2)
420               ENDIF
421              ENDDO
422             ENDDO
423           ENDIF
424    
425  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
426  CADJ STORE tsf  = comlev1_thsice_3, key=ikey_3  CADJ STORE devdt(i,j) = comlev1_bibj,key=iicekey,byte=isbyte
427    CADJ STORE dflxdt(i,j) = comlev1_bibj,key=iicekey,byte=isbyte
428    CADJ STORE flxtexsw(i,j) = comlev1_bibj,key=iicekey,byte=isbyte
429    CADJ STORE iceflag(i,j) = comlev1_bibj,key=iicekey,byte=isbyte
430    CADJ STORE tsf(i,j) = comlev1_bibj,key=iicekey,byte=isbyte
431  #endif  #endif
432           a1 = a10 - k12*df0dT / (k12-df0dT)  
433           b1 = b10 - k12*(flx0-df0dT*Tsf) / (k12-df0dT)  C--   Compute new top layer and surface temperatures.
434           Tice(1) = -(b1 + SQRT(b1*b1-4. _d 0*a1*c1))/(2. _d 0*a1)  C     If Tsfc is computed to be > 0 C, fix Tsfc = 0 and recompute T1
435           dTsf = (flx0 + k12*(Tice(1)-Tsf)) / (k12-df0dT)  C     with different coefficients.
436           Tsf = Tsf + dTsf         DO j = jMin, jMax
437           IF (Tsf .GT. 0. _d 0) THEN          DO i = iMin, iMax
438             IF ( iceFlag(i,j) ) THEN
439              flxNet = sHeat(i,j) + flxTexSW(i,j)
440  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
441              IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)            IF ( dBug(i,j,bi,bj) ) WRITE(6,1020)
442       &        'ThSI_SOLVE4T: k,ts,t1,dTs=', k,Tsf,Tice(1),dTsf       &     'ThSI_SOLVE4T: flxNet,dFlxdT,k12,D=',
443         &      flxNet, dFlxdT(i,j), k12(i,j), k12(i,j)-dFlxdT(i,j)
444  #endif  #endif
445              a1 = a10 + k12  
446              b1 = b10          ! note b1 = b10 - k12*Tf0            a1 = a10(i,j) - k12(i,j)*dFlxdT(i,j) / (k12(i,j)-dFlxdT(i,j))
447              Tice(1) = (-b1 - SQRT(b1*b1-4. _d 0*a1*c1))/(2. _d 0*a1)            b1 = b10(i,j) - k12(i,j)*(flxNet-dFlxdT(i,j)*Tsf(i,j))
448              Tsf = 0. _d 0       &                  /(k12(i,j)-dFlxdT(i,j))
449              c1 = c10(i,j)
450              tIc1(i,j)  = -(b1 + SQRT(b1*b1-4. _d 0*a1*c1))/(2. _d 0*a1)
451              dTsrf(i,j) = (flxNet + k12(i,j)*(tIc1(i,j)-Tsf(i,j)))
452         &                  /(k12(i,j)-dFlxdT(i,j))
453              Tsf(i,j) = Tsf(i,j) + dTsrf(i,j)
454              IF ( Tsf(i,j) .GT. 0. _d 0 ) THEN
455    #ifdef ALLOW_DBUG_THSICE
456               IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)
457         &     'ThSI_SOLVE4T: k,ts,t1,dTs=',k,Tsf(i,j),tIc1(i,j),dTsrf(i,j)
458    #endif
459               a1 = a10(i,j) + k12(i,j)
460    C      note: b1 = b10 - k12*Tf0
461               b1 = b10(i,j)
462               tIc1(i,j) = (-b1 - SQRT(b1*b1-4. _d 0*a1*c1))/(2. _d 0*a1)
463               Tsf(i,j) = 0. _d 0
464             IF ( useBlkFlx ) THEN             IF ( useBlkFlx ) THEN
465              IF (hs.GT.3. _d -1) THEN              flxTexSW(i,j) = flx0exSW(i,j)
466                   iceornot=2              evapT(i,j) = evap0(i,j)
467              ELSE              dTsrf(i,j) = 0. _d 0
                  iceornot=1  
             ENDIF  
             IF ( useBulkForce ) THEN  
              CALL THSICE_GET_BULKF(  
      I                            iceornot, Tsf,  
      O                            flxExceptSw, df0dT, evap, dEvdT,  
      I                            i,j,bi,bj,myThid )  
             ELSEIF ( useEXF ) THEN  
              CALL THSICE_GET_EXF  (  
      I                            iceornot, Tsf,  
      O                            flxExceptSw, df0dT, evap, dEvdT,  
      I                            i,j,bi,bj,myThid )  
             ENDIF  
             dTsf = 0. _d 0  
468             ELSE             ELSE
469              flxExceptSw = flxExSW(i,j,0)              flxTexSW(i,j) = flxExSW(i,j,0)
470              dTsf = 1000.              dTsrf(i,j) = 1000.
471              df0dT = 0.              dFlxdT(i,j) = 0.
472             ENDIF             ENDIF
473             flx0 = fswdn + flxExceptSw            ENDIF
          ENDIF  
474    
475  C Check for convergence.  If no convergence, then repeat.  C--   Check for convergence.  If no convergence, then repeat.
476  C            iceFlag(i,j) = ABS(dTsrf(i,j)).GE.Terrmax
477  C Convergence test: Make sure Tsfc has converged, within prescribed error.            iterate4Tsf = iterate4Tsf .OR. iceFlag(i,j)
478  C (Energy conservation is guaranteed within machine roundoff, even  
479  C if Tsfc has not converged.)  C     Convergence test: Make sure Tsfc has converged, within prescribed error.
480  C If no convergence, then repeat.  C     (Energy conservation is guaranteed within machine roundoff, even
481    C      if Tsfc has not converged.)
482    C     If no convergence, then repeat.
483    
484  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
485           IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)            IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)
486       &     'ThSI_SOLVE4T: k,ts,t1,dTs=', k,Tsf,Tice(1),dTsf       &    'ThSI_SOLVE4T: k,ts,t1,dTs=', k,Tsf(i,j),tIc1(i,j),dTsrf(i,j)
487              IF ( useBlkFlx .AND. k.EQ.nitMaxTsf .AND. iceFlag(i,j) ) THEN
488               WRITE (6,'(A,4I4,I12,F15.9)')
489         &       ' BB: not converge: i,j,it,hi=',i,j,bi,bj,myIter,hIce(i,j)
490               WRITE (6,*) 'BB: not converge: Tsf, dTsf=',
491         &          Tsf(i,j), dTsrf(i,j)
492               WRITE (6,*) 'BB: not converge: flxNet,dFlxT=',
493         &          flxNet, dFlxdT(i,j)
494               IF (Tsf(i,j).LT.-70. _d 0) STOP
495              ENDIF
496  #endif  #endif
          IF ( useBlkFlx .AND. k.EQ.nitMaxTsf  
      &                  .AND. ABS(dTsf).GE.Terrmax ) THEN  
             WRITE (6,'(A,4I4,I12,F15.9)')  
      &                 ' BB: not converge: i,j,it,hi=',i,j,bi,bj,  
      &                                                    myIter,hi  
             WRITE (6,*) 'BB: not converge: Tsf, dTsf=', Tsf,dTsf  
             WRITE (6,*) 'BB: not converge: flx0,dfdT=',flx0,df0dT  
             IF (Tsf.LT.-70. _d 0) STOP  
          ENDIF  
497    
498             ENDIF
499            ENDDO
500           ENDDO
501    #ifndef ALLOW_AUTODIFF_TAMC
502         ENDIF         ENDIF
503    #endif
504        ENDDO        ENDDO
505  C ------ end iteration ------------  C ------ end iteration ------------
506    
507  C Compute new bottom layer temperature.  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|
508    
509  #ifdef ALLOW_AUTODIFF_TAMC        DO j = jMin, jMax
510  CADJ STORE  Tice(:)      = comlev1_thsice_1, key=ikey_1         DO i = iMin, iMax
511  CADJ STORE  df0dt        = comlev1_thsice_1, key=ikey_1          IF ( iceMask(i,j).GT.0. _d 0 ) THEN
512  #endif  
513        Tice(2) = (2. _d 0*dt*k32*(Tice(1)+2. _d 0*Tf)  C--   Compute new bottom layer temperature.
514       &        + rhoi*cpice *hi*Tice(2))            k32       = 2. _d 0*kIce  / hIce(i,j)
515       &         /(6. _d 0*dt*k32 + rhoi*cpice *hi)            tIc2(i,j) = ( 2. _d 0*dt*k32*(tIc1(i,j)+2. _d 0*tFrz(i,j))
516         &                 + rhoi*cpIce*hIce(i,j)*tIc2(i,j))
517         &               /(6. _d 0*dt*k32 + rhoi*cpIce*hIce(i,j))
518  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
519        IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)            IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)
520       &  'ThSI_SOLVE4T: k, Ts, Tice=',k,Tsf,Tice       &  'ThSI_SOLVE4T: k, Ts, Tice=',k,Tsf(i,j),tIc1(i,j),tIc2(i,j)
521              netSW   = flxAtm(i,j)
522  #endif  #endif
523    
524  C Compute final flux values at surfaces.  C--   Compute final flux values at surfaces.
525              tSrf(i,j)  = Tsf(i,j)
526        fct    = k12*(Tsf-Tice(1))            fct    = k12(i,j)*(Tsf(i,j)-tIc1(i,j))
527        flxCnB(i,j) = 4. _d 0*kice *(Tice(2)-Tf)/hi            flxCnB(i,j) = 4. _d 0*kIce *(tIc2(i,j)-tFrz(i,j))/hIce(i,j)
528        flx0   = flx0 + df0dT*dTsf            flxNet = sHeat(i,j) + flxTexSW(i,j)
529        IF ( useBlkFlx ) THEN            flxNet = flxNet + dFlxdT(i,j)*dTsrf(i,j)
530  C--   needs to update also Evap (Tsf changes) since Latent heat has been updated            IF ( useBlkFlx ) THEN
531          evpAtm(i,j) = evap + dEvdT*dTsf  C-    needs to update also Evap (Tsf changes) since Latent heat has been updated
532        ELSE              evpAtm(i,j) = evapT(i,j) + dEvdT(i,j)*dTsrf(i,j)
533              ELSE
534  C- WARNING: Evap & +Evap*Lfresh are missing ! (but only affects Diagnostics)  C- WARNING: Evap & +Evap*Lfresh are missing ! (but only affects Diagnostics)
535          evpAtm(i,j) = 0.              evpAtm(i,j) = 0.
536        ENDIF            ENDIF
537  C-    energy flux to Atmos: use net short-wave flux at surf. and  C-    Update energy flux to Atmos with other than SW contributions;
538  C     use latent heat = Lvap (energy=0 for liq. water at 0.oC)  C     use latent heat = Lvap (energy=0 for liq. water at 0.oC)
539        flxAtm(i,j) = netSW + flxExceptSw            flxAtm(i,j) = flxAtm(i,j) + flxTexSW(i,j)
540       &                    + df0dT*dTsf + evpAtm(i,j)*Lfresh       &                + dFlxdT(i,j)*dTsrf(i,j) + evpAtm(i,j)*Lfresh
541  C-    excess of energy @ surface (used for surface melting):  C-    excess of energy @ surface (used for surface melting):
542        sHeat(i,j) = flx0 - fct            sHeat(i,j) = flxNet - fct
   
 C-    SW flux at sea-ice base left to the ocean  
       flxSW(i,j) = fswocn  
543    
544  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
545        IF ( dBug(i,j,bi,bj) ) WRITE(6,1020)            IF ( dBug(i,j,bi,bj) ) WRITE(6,1020)
546       &  'ThSI_SOLVE4T: flx0,fct,Dif,flxCnB=',       &     'ThSI_SOLVE4T: flxNet,fct,Dif,flxCnB=',
547       &                 flx0,fct,flx0-fct,flxCnB(i,j)       &                    flxNet,fct,flxNet-fct,flxCnB(i,j)
548  #endif  #endif
549    
550  C Compute new enthalpy for each layer.  C--   Compute new enthalpy for each layer.
551              qIc1(i,j) = -cpWater*Tmlt1 + cpIce *(Tmlt1-tIc1(i,j))
552         &                + Lfresh*(1. _d 0-Tmlt1/tIc1(i,j))
553              qIc2(i,j) = -cpIce *tIc2(i,j) + Lfresh
554    
       qicen(1) = -cpwater*Tmlt1 + cpice *(Tmlt1-Tice(1))  
      &            + Lfresh*(1. _d 0-Tmlt1/Tice(1))  
       qicen(2) = -cpice *Tice(2) + Lfresh  
   
 C Make sure internal ice temperatures do not exceed Tmlt.  
 C (This should not happen for reasonable values of i0.)  
   
       IF (Tice(1) .GE. Tmlt1) THEN  
         WRITE (6,'(A,2I4,2I3,1P2E14.6)')  
      &   ' BBerr - Bug: IceT(1) > Tmlt',i,j,bi,bj,Tice(1),Tmlt1  
       ENDIF  
       IF (Tice(2) .GE. 0. _d 0) THEN  
        WRITE (6,'(A,2I4,2I3,1P2E14.6)')  
      &   ' BBerr - Bug: IceT(2) > 0',i,j,bi,bj,Tice(2)  
       ENDIF  
   
 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  
 C--    Update Sea-Ice state :  
           tSrf(i,j)  = Tsf  
           tIc1(i,j)  = Tice(1)  
           tic2(i,j)  = Tice(2)  
           qIc1(i,j)  = qicen(1)  
           qIc2(i,j)  = qicen(2)  
 c         dTsrf(i,j) = dTsf  
           IF ( .NOT.useBlkFlx ) dTsrf(i,j) = dTsf  
 c         sHeat(i,j) = sHeating  
 c         flxCnB(i,j)= flxCnB  
 c         flxAtm(i,j)= flxAtm  
 c         evpAtm(i,j)= evpAtm  
555  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
556    C--   Make sure internal ice temperatures do not exceed Tmlt.
557    C     (This should not happen for reasonable values of i0swFrac)
558              IF (tIc1(i,j) .GE. Tmlt1) THEN
559               WRITE (6,'(A,2I4,2I3,1P2E14.6)')
560         &     ' BBerr - Bug: IceT(1) > Tmlt',i,j,bi,bj,tIc1(i,j),Tmlt1
561              ENDIF
562              IF (tIc2(i,j) .GE. 0. _d 0) THEN
563               WRITE (6,'(A,2I4,2I3,1P2E14.6)')
564         &     ' BBerr - Bug: IceT(2) > 0',i,j,bi,bj,tIc2(i,j)
565              ENDIF
566    
567            IF ( dBug(i,j,bi,bj) ) THEN            IF ( dBug(i,j,bi,bj) ) THEN
568             WRITE(6,1020) 'ThSI_SOLV_4T: Tsf, Tice(1,2), dTsurf=',             WRITE(6,1020) 'ThSI_SOLV_4T: Tsf, Tice(1,2), dTsurf=',
569       &                              Tsf, Tice, dTsf       &           Tsf(i,j), tIc1(i,j), tIc2(i,j), dTsrf(i,j)
570             WRITE(6,1020) 'ThSI_SOLV_4T: sHeat, flxCndBt, Qice =',             WRITE(6,1020) 'ThSI_SOLV_4T: sHeat, flxCndBt, Qice =',
571       &           sHeat(i,j), flxCnB(i,j), qicen       &           sHeat(i,j), flxCnB(i,j), qIc1(i,j), qIc2(i,j)
572             WRITE(6,1020) 'ThSI_SOLV_4T: flxA, evpA, fxSW_bf,af=',             WRITE(6,1020) 'ThSI_SOLV_4T: flxA, evpA, fxSW_bf,af=',
573       &           flxAtm(i,j), evpAtm(i,j), netSW, flxSW(i,j)       &           flxAtm(i,j), evpAtm(i,j), netSW, flxSW(i,j)
574            ENDIF            ENDIF
575  #endif  #endif
576    
577          ELSE          ELSE
578            IF ( .NOT.useBlkFlx ) dTsrf(i,j) = 0. _d 0  C--     ice-free grid point:
579    c         tIc1  (i,j) = 0. _d 0
580    c         tIc2  (i,j) = 0. _d 0
581              dTsrf (i,j) = 0. _d 0
582    c         sHeat (i,j) = 0. _d 0
583    c         flxCnB(i,j) = 0. _d 0
584    c         flxAtm(i,j) = 0. _d 0
585    c         evpAtm(i,j) = 0. _d 0
586    
587          ENDIF          ENDIF
588         ENDDO         ENDDO
589        ENDDO        ENDDO
590  #endif  /* ALLOW_THSICE */  #endif  /* ALLOW_THSICE */
591    
592  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|
593    
594        RETURN        RETURN
595        END        END
Legend:
Removed from v.1.15  
changed lines
  Added in v.1.24
  ViewVC Help
Powered by ViewVC 1.1.22