/[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.17 by jmc, Tue May 1 02:26:50 2007 UTC revision 1.27 by jmc, Fri Oct 22 13:43:35 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,
# 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        _RL  recip_dhSnowLin  C     a10,b10,c10  :: coefficients in quadratic eqn for T1
140        INTEGER iceornot  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, k, iterMax
146          INTEGER ii, jj, icount, stdUnit
147          CHARACTER*(MAX_LEN_MBUF) msgBuf
148          _RL  frsnow
149          _RL  fswpen
150          _RL  fswdn
151          _RL  fswint
152          _RL  fswocn
153          _RL  Tsf     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
154          _RL  flx0exSW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
155          _RL  flxTexSW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
156          _RL  dFlxdT  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
157          _RL  evap0   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
158          _RL  evapT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
159          _RL  dEvdT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
160          _RL  flxNet
161          _RL  fct
162          _RL  k12     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
163          _RL  k32
164          _RL  a10     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
165          _RL  b10     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
166          _RL  c10     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
167          _RL  a1, b1, c1
168          _RL  dt
169          _RL  recip_dhSnowLin
170    #ifdef ALLOW_DBUG_THSICE
171          _RL  netSW
172    #endif
173    
174  C-    define grid-point location where to print debugging values  C-    Define grid-point location where to print debugging values
175  #include "THSICE_DEBUG.h"  #include "THSICE_DEBUG.h"
176    
177   1010 FORMAT(A,I3,3F11.6)   1010 FORMAT(A,I3,3F11.6)
178   1020 FORMAT(A,1P4E14.6)   1020 FORMAT(A,1P4E14.6)
179    
180  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|
181    
182  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
183        act1 = bi - myBxLo(myThid)        act1 = bi - myBxLo(myThid)
# Line 175  C---+----1----+----2----+----3----+----4 Line 187  C---+----1----+----2----+----3----+----4
187        act3 = myThid - 1        act3 = myThid - 1
188        max3 = nTx*nTy        max3 = nTx*nTy
189        act4 = ikey_dynamics - 1        act4 = ikey_dynamics - 1
190          iicekey = (act1 + 1) + act2*max1
191         &                     + act3*max1*max2
192         &                     + act4*max1*max2*max3
193  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
194    
195    #ifdef ALLOW_AUTODIFF_TAMC
196    CADJ STORE flxsw(:,:) = comlev1_bibj,key=iicekey,byte=isbyte
197          DO j = 1-OLy, sNy+OLy
198           DO i = 1-OLx, sNx+OLx
199            tIc1(i,j) = 0. _d 0
200            tIc2(i,j) = 0. _d 0
201           ENDDO
202          ENDDO
203    #endif
204    
205          stdUnit = standardMessageUnit
206        useBlkFlx = useEXF .OR. useBulkForce        useBlkFlx = useEXF .OR. useBulkForce
207          dt  = thSIce_dtTemp
208        IF ( dhSnowLin.GT.0. _d 0 ) THEN        IF ( dhSnowLin.GT.0. _d 0 ) THEN
209          recip_dhSnowLin = 1. _d 0 / dhSnowLin          recip_dhSnowLin = 1. _d 0 / dhSnowLin
210        ELSE        ELSE
211          recip_dhSnowLin = 0. _d 0          recip_dhSnowLin = 0. _d 0
212        ENDIF        ENDIF
213    
214        dt  = thSIce_dtTemp        iterate4Tsf = .FALSE.
215          icount = 0
216    
217        DO j = jMin, jMax        DO j = jMin, jMax
218         DO i = iMin, iMax         DO i = iMin, iMax
219  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
220            ikey_1 = i  c         ikey_1 = i
221       &         + sNx*(j-1)  c    &         + sNx*(j-1)
222       &         + sNx*sNy*act1  c    &         + sNx*sNy*act1
223       &         + sNx*sNy*max1*act2  c    &         + sNx*sNy*max1*act2
224       &         + sNx*sNy*max1*max2*act3  c    &         + sNx*sNy*max1*max2*act3
225       &         + sNx*sNy*max1*max2*max3*act4  c    &         + sNx*sNy*max1*max2*max3*act4
226  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
 C--  
227  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
228  CADJ STORE  devdt        = comlev1_thsice_1, key=ikey_1  cCADJ STORE  devdt(i,j)        = comlev1_thsice_1, key=ikey_1
229  CADJ STORE  df0dt        = comlev1_thsice_1, key=ikey_1  cCADJ STORE  dFlxdT(i,j)        = comlev1_thsice_1, key=ikey_1
230  CADJ STORE  flxexceptsw  = comlev1_thsice_1, key=ikey_1  cCADJ STORE  flxsw(i,j)   = comlev1_thsice_1, key=ikey_1
231  CADJ STORE  flxsw(i,j)   = comlev1_thsice_1, key=ikey_1  cCADJ STORE  qic1(i,j)    = comlev1_thsice_1, key=ikey_1
232  CADJ STORE  qic1(i,j)    = comlev1_thsice_1, key=ikey_1  cCADJ STORE  qic2(i,j)    = comlev1_thsice_1, key=ikey_1
233  CADJ STORE  qic2(i,j)    = comlev1_thsice_1, key=ikey_1  cCADJ STORE  tsrf(i,j)    = comlev1_thsice_1, key=ikey_1
 CADJ STORE  tsrf(i,j)    = comlev1_thsice_1, key=ikey_1  
234  #endif  #endif
235          IF ( iceMask(i,j).GT.0. _d 0) THEN          IF ( iceMask(i,j).GT.0. _d 0) THEN
236            hi      = hIce(i,j)            iterate4Tsf  = .TRUE.
237            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-|--+----|  
238  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
239            IF ( dBug(i,j,bi,bj) ) WRITE(6,'(A,2I4,2I2)')            IF ( dBug(i,j,bi,bj) ) WRITE(stdUnit,'(A,2I4,2I2)')
240       &         'ThSI_SOLVE4T: i,j=',i,j,bi,bj       &         'ThSI_SOLVE4T: i,j=',i,j,bi,bj
241  #endif  #endif
242            IF ( hi.LT.hIceMin ) THEN            IF ( hIce(i,j).LT.hIceMin ) THEN
243  C If hi < hIceMin, melt the ice.  C     if hi < hIceMin, melt the ice.
244               STOP 'THSICE_SOLVE4TEMP: should not enter if hi<hIceMin'  C     keep the position of this problem but do the stop later
245               ii = i
246               jj = j
247               icount = icount + 1
248            ENDIF            ENDIF
249    
250  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C--   Fractional snow cover:
251    C     assume a linear distribution of snow thickness, with dhSnowLin slope,
252    C      from hs-dhSnowLin to hs+dhSnowLin if full ice & snow cover.
253    C     frsnow = fraction of snow over the ice-covered part of the grid cell
254              IF ( hSnow(i,j) .GT. iceMask(i,j)*dhSnowLin ) THEN
255               frsnow = 1. _d 0
256              ELSE
257               frsnow = hSnow(i,j)*recip_dhSnowLin/iceMask(i,j)
258               IF ( frsnow.GT.0. _d 0 ) frsnow = SQRT(frsnow)
259              ENDIF
260    
261  C fractional snow cover  C--   Compute SW flux absorbed at surface and penetrating to layer 1.
262  C  assume a linear distribution of snow thickness, with dhSnowLin slope,            fswpen = flxSW(i,j) * (1. _d 0 - frsnow) * i0swFrac
263  C   from hs-dhSnowLin to hs+dhSnowLin if full ice & snow cover.            fswocn = fswpen * exp(-ksolar*hIce(i,j))
264  C frsnow = fraction of snow over the ice-covered part of the grid cell            fswint = fswpen - fswocn
265        IF ( hs .GT. iceMask(i,j)*dhSnowLin ) THEN            fswdn  = flxSW(i,j) - fswpen
266          frsnow = 1. _d 0  
267        ELSE  C     Initialise Atmospheric surf. heat flux with net SW flux at the surface
268          frsnow = hs*recip_dhSnowLin/iceMask(i,j)            flxAtm(i,j) = flxSW(i,j)
269          IF ( frsnow.GT.0. _d 0 ) frsnow = SQRT(frsnow)  C     SW flux at sea-ice base left to the ocean
270        ENDIF            flxSW(i,j) = fswocn
271    C     Initialise surface Heating with SW contribution
272              sHeat(i,j) = fswdn
273    
274    C--   Compute conductivity terms at layer interfaces.
275              k12(i,j) = 4. _d 0*kIce*kSnow
276         &             / (kSnow*hIce(i,j) + 4. _d 0*kIce*hSnow(i,j))
277              k32      = 2. _d 0*kIce  / hIce(i,j)
278    
279    C--   Compute ice temperatures
280              a1 = cpIce
281              b1 = qIc1(i,j) + (cpWater-cpIce )*Tmlt1 - Lfresh
282              c1 = Lfresh * Tmlt1
283              tIc1(i,j) = 0.5 _d 0 *(-b1 - SQRT(b1*b1-4. _d 0*a1*c1))/a1
284              tIc2(i,j) = (Lfresh-qIc2(i,j)) / cpIce
285    
 C Compute SW flux absorbed at surface and penetrating to layer 1.  
       fswpen  = netSW * (1. _d 0 - frsnow) * i0swFrac  
       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  
286  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
287        IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)            IF (tIc1(i,j).GT.0. _d 0 ) THEN
288       &  'ThSI_SOLVE4T: k, Ts, Tice=',0,Tsf,Tice             WRITE(stdUnit,'(A,I12,1PE14.6)')
289         &       ' BBerr: Tice(1) > 0 ; it=', myIter, qIc1(i,j)
290               WRITE(stdUnit,'(A,4I5,2F11.4)')
291         &      ' BBerr: i,j,bi,bj,Tice = ',i,j,bi,bj,tIc1(i,j),tIc2(i,j)
292              ENDIF
293              IF ( tIc2(i,j).GT.0. _d 0) THEN
294               WRITE(stdUnit,'(A,I12,1PE14.6)')
295         &       ' BBerr: Tice(2) > 0 ; it=', myIter, qIc2(i,j)
296               WRITE(stdUnit,'(A,4I5,2F11.4)')
297         &      ' BBerr: i,j,bi,bj,Tice = ',i,j,bi,bj,tIc1(i,j),tIc2(i,j)
298              ENDIF
299              IF ( dBug(i,j,bi,bj) ) WRITE(stdUnit,1010)
300         &      'ThSI_SOLVE4T: k, Ts, Tice=',0,tSrf(i,j),tIc1(i,j),tIc2(i,j)
301  #endif  #endif
302    
303  C Compute coefficients used in quadratic formula.  C--   Compute coefficients used in quadratic formula.
304    
305        a10 = rhoi*cpIce *hi/(2. _d 0*dt) +            a10(i,j) = rhoi*cpIce *hIce(i,j)/(2. _d 0*dt) +
306       &      k32 * (4. _d 0*dt*k32 + rhoi*cpIce *hi)       &          k32*( 4. _d 0*dt*k32 + rhoi*cpIce *hIce(i,j) )
307       &       / (6. _d 0*dt*k32 + rhoi*cpIce *hi)       &           / ( 6. _d 0*dt*k32 + rhoi*cpIce *hIce(i,j) )
308        b10 = -hi*            b10(i,j) = -hIce(i,j)*
309       &      (rhoi*cpIce*Tice(1)+rhoi*Lfresh*Tmlt1/Tice(1))       &          ( rhoi*cpIce*tIc1(i,j) + rhoi*Lfresh*Tmlt1/tIc1(i,j) )
310       &       /(2. _d 0*dt)       &           /(2. _d 0*dt)
311       &      - k32 * (4. _d 0*dt*k32*Tf+rhoi*cpIce *hi*Tice(2))       &        - k32*( 4. _d 0*dt*k32*tFrz(i,j)
312       &       / (6. _d 0*dt*k32 + rhoi*cpIce *hi) - fswint       &               +rhoi*cpIce*hIce(i,j)*tIc2(i,j) )
313        c1 = rhoi*Lfresh*hi*Tmlt1 / (2. _d 0*dt)       &           / ( 6. _d 0*dt*k32 + rhoi*cpIce *hIce(i,j) )
314         &        - fswint
315  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.  
316    
317        IF ( useBlkFlx ) THEN          ELSE
318          iterMax = nitMaxTsf            iceFlag(i,j) = .FALSE.
319        ELSE          ENDIF
320          iterMax = 1         ENDDO
321          ENDDO
322          IF ( icount .gt. 0 ) THEN
323           WRITE(stdUnit,'(A,I5,A)')
324         &      'THSICE_SOLVE4TEMP: there are ',icount,
325         &      ' case(s) where hIce<hIceMin;'
326           WRITE(stdUnit,'(A,I3,A1,I3,A)')
327         &      'THSICE_SOLVE4TEMP: the last one was at (',ii,',',jj,
328         &      ') with hIce = ', hIce(ii,jj)
329           WRITE( msgBuf, '(A)')
330         &      'THSICE_SOLVE4TEMP: should not enter if hIce<hIceMin'
331           CALL PRINT_ERROR( msgBuf , myThid )
332           STOP 'ABNORMAL END: S/R THSICE_SOLVE4TEMP'
333        ENDIF        ENDIF
       dTsf = Terrmax  
334    
335  C ----- begin iteration  -----  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|
       DO k = 1,iterMax  
336    
337  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
338         ikey_3 = (ikey_1-1)*MaxTsf + k  CADJ STORE devdt(:,:)    = comlev1_bibj,key=iicekey,byte=isbyte
339    CADJ STORE tsf(:,:)      = comlev1_bibj,key=iicekey,byte=isbyte
340  #endif  #endif
341    
342    C--   Get surface fluxes over melting surface
343  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
344  CADJ STORE tsf          = comlev1_thsice_3, key=ikey_3        IF ( useBlkFlx ) THEN
345  CADJ STORE dtsf         = comlev1_thsice_3, key=ikey_3  #else
346  CADJ STORE df0dt        = comlev1_thsice_3, key=ikey_3        IF ( useBlkFlx .AND. iterate4Tsf  ) THEN
 CADJ STORE flxexceptsw  = comlev1_thsice_3, key=ikey_3  
347  #endif  #endif
348         IF ( ABS(dTsf).GE.Terrmax ) THEN          DO j = jMin, jMax
349             DO i = iMin, iMax
350               Tsf(i,j) = 0.
351             ENDDO
352            ENDDO
353            IF ( useEXF ) THEN
354               CALL THSICE_GET_EXF(   bi, bj,
355         I                            iMin, iMax, jMin, jMax,
356         I                            iceFlag, hSnow, Tsf,
357         O                            flx0exSW, dFlxdT, evap0, dEvdT,
358         I                            myTime, myIter, myThid )
359    C-    could add this "ifdef" to hide THSICE_GET_BULKF from TAF
360    #ifdef ALLOW_BULK_FORCE
361            ELSEIF ( useBulkForce ) THEN
362               CALL THSICE_GET_BULKF( bi, bj,
363         I                            iMin, iMax, jMin, jMax,
364         I                            iceFlag, hSnow, Tsf,
365         O                            flx0exSW, dFlxdT, evap0, dEvdT,
366         I                            myTime, myIter, myThid )
367    #endif /* ALLOW_BULK_FORCE */
368            ENDIF
369          ENDIF
370    
371  C Save temperatures at start of iteration.  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|
372  c        Tsf_start = Tsf  C---  Compute new surface and internal temperatures; iterate until
373    C     Tsfc converges.
374          DO j = jMin, jMax
375            DO i = iMin, iMax
376              Tsf(i,j)  = tSrf(i,j)
377              dTsrf(i,j) = Terrmax
378            ENDDO
379          ENDDO
380          IF ( useBlkFlx ) THEN
381            iterMax = nitMaxTsf
382          ELSE
383            iterMax = 1
384          ENDIF
385    
         IF ( useBlkFlx ) THEN  
386  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
387  CADJ STORE tsf          = comlev1_thsice_3, key=ikey_3  CADJ STORE devdt(:,:)    = comlev1_bibj,key=iicekey,byte=isbyte
388  #endif  CADJ STORE dflxdt(:,:)   = comlev1_bibj,key=iicekey,byte=isbyte
389  C Compute top surface flux.  CADJ STORE flx0exsw(:,:) = comlev1_bibj,key=iicekey,byte=isbyte
390           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  
391  #endif  #endif
392    
393  C Compute new top layer and surface temperatures.  C ----- begin iteration  -----
394  C If Tsfc is computed to be > 0 C, fix Tsfc = 0 and recompute T1        DO k = 1,iterMax
395  C with different coefficients.  #ifndef ALLOW_AUTODIFF_TAMC
396           IF ( iterate4Tsf ) THEN
397           iterate4Tsf = .FALSE.
398    #endif
399    
400           IF ( useBlkFlx ) THEN
401    C--   Compute top surface flux.
402             IF ( useEXF ) THEN
403               CALL THSICE_GET_EXF(   bi, bj,
404         I                            iMin, iMax, jMin, jMax,
405         I                            iceFlag, hSnow, Tsf,
406         O                            flxTexSW, dFlxdT, evapT, dEvdT,
407         I                            myTime, myIter, myThid )
408    C-    could add this "ifdef" to hide THSICE_GET_BULKF from TAF
409    #ifdef ALLOW_BULK_FORCE
410             ELSEIF ( useBulkForce ) THEN
411               CALL THSICE_GET_BULKF( bi, bj,
412         I                            iMin, iMax, jMin, jMax,
413         I                            iceFlag, hSnow, Tsf,
414         O                            flxTexSW, dFlxdT, evapT, dEvdT,
415         I                            myTime, myIter, myThid )
416    #endif /* ALLOW_BULK_FORCE */
417             ENDIF
418           ELSE
419             DO j = jMin, jMax
420              DO i = iMin, iMax
421               IF ( iceFlag(i,j) ) THEN
422                 flxTexSW(i,j) = flxExSW(i,j,1)
423                 dFlxdT(i,j) = flxExSW(i,j,2)
424               ENDIF
425              ENDDO
426             ENDDO
427           ENDIF
428    
429  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
430  CADJ STORE tsf  = comlev1_thsice_3, key=ikey_3  CADJ STORE devdt(:,:)    = comlev1_bibj,key=iicekey,byte=isbyte
431    CADJ STORE dflxdt(:,:)   = comlev1_bibj,key=iicekey,byte=isbyte
432    CADJ STORE flxtexsw(:,:) = comlev1_bibj,key=iicekey,byte=isbyte
433    CADJ STORE iceflag(:,:)  = comlev1_bibj,key=iicekey,byte=isbyte
434    CADJ STORE tsf(:,:)      = comlev1_bibj,key=iicekey,byte=isbyte
435  #endif  #endif
436           a1 = a10 - k12*df0dT / (k12-df0dT)  
437           b1 = b10 - k12*(flx0-df0dT*Tsf) / (k12-df0dT)  C--   Compute new top layer and surface temperatures.
438           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
439           dTsf = (flx0 + k12*(Tice(1)-Tsf)) / (k12-df0dT)  C     with different coefficients.
440           Tsf = Tsf + dTsf         DO j = jMin, jMax
441           IF (Tsf .GT. 0. _d 0) THEN          DO i = iMin, iMax
442             IF ( iceFlag(i,j) ) THEN
443              flxNet = sHeat(i,j) + flxTexSW(i,j)
444  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
445              IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)            IF ( dBug(i,j,bi,bj) ) WRITE(stdUnit,1020)
446       &        'ThSI_SOLVE4T: k,ts,t1,dTs=', k,Tsf,Tice(1),dTsf       &     'ThSI_SOLVE4T: flxNet,dFlxdT,k12,D=',
447         &      flxNet, dFlxdT(i,j), k12(i,j), k12(i,j)-dFlxdT(i,j)
448  #endif  #endif
449              a1 = a10 + k12  
450              b1 = b10          ! note b1 = b10 - k12*Tf0            a1 = a10(i,j) - k12(i,j)*dFlxdT(i,j) / (k12(i,j)-dFlxdT(i,j))
451              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))
452              Tsf = 0. _d 0       &                  /(k12(i,j)-dFlxdT(i,j))
453              c1 = c10(i,j)
454              tIc1(i,j)  = -(b1 + SQRT(b1*b1-4. _d 0*a1*c1))/(2. _d 0*a1)
455              dTsrf(i,j) = (flxNet + k12(i,j)*(tIc1(i,j)-Tsf(i,j)))
456         &                  /(k12(i,j)-dFlxdT(i,j))
457              Tsf(i,j) = Tsf(i,j) + dTsrf(i,j)
458              IF ( Tsf(i,j) .GT. 0. _d 0 ) THEN
459    #ifdef ALLOW_DBUG_THSICE
460               IF ( dBug(i,j,bi,bj) ) WRITE(stdUnit,1010)
461         &     'ThSI_SOLVE4T: k,ts,t1,dTs=',k,Tsf(i,j),tIc1(i,j),dTsrf(i,j)
462    #endif
463               a1 = a10(i,j) + k12(i,j)
464    C      note: b1 = b10 - k12*Tf0
465               b1 = b10(i,j)
466               tIc1(i,j) = (-b1 - SQRT(b1*b1-4. _d 0*a1*c1))/(2. _d 0*a1)
467               Tsf(i,j) = 0. _d 0
468             IF ( useBlkFlx ) THEN             IF ( useBlkFlx ) THEN
469              IF (hs.GT.3. _d -1) THEN              flxTexSW(i,j) = flx0exSW(i,j)
470                   iceornot=2              evapT(i,j) = evap0(i,j)
471              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  
472             ELSE             ELSE
473              flxExceptSw = flxExSW(i,j,0)              flxTexSW(i,j) = flxExSW(i,j,0)
474              dTsf = 1000.              dTsrf(i,j) = 1000.
475              df0dT = 0.              dFlxdT(i,j) = 0.
476             ENDIF             ENDIF
477             flx0 = fswdn + flxExceptSw            ENDIF
          ENDIF  
478    
479  C Check for convergence.  If no convergence, then repeat.  C--   Check for convergence.  If no convergence, then repeat.
480  C            iceFlag(i,j) = ABS(dTsrf(i,j)).GE.Terrmax
481  C Convergence test: Make sure Tsfc has converged, within prescribed error.            iterate4Tsf = iterate4Tsf .OR. iceFlag(i,j)
482  C (Energy conservation is guaranteed within machine roundoff, even  
483  C if Tsfc has not converged.)  C     Convergence test: Make sure Tsfc has converged, within prescribed error.
484  C If no convergence, then repeat.  C     (Energy conservation is guaranteed within machine roundoff, even
485    C      if Tsfc has not converged.)
486    C     If no convergence, then repeat.
487    
488  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
489           IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)            IF ( dBug(i,j,bi,bj) ) WRITE(stdUnit,1010)
490       &     '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)
491              IF ( useBlkFlx .AND. k.EQ.nitMaxTsf .AND. iceFlag(i,j) ) THEN
492               WRITE(stdUnit,'(A,4I4,I12,F15.9)')
493         &       ' BB: not converge: i,j,it,hi=',i,j,bi,bj,myIter,hIce(i,j)
494               WRITE(stdUnit,*)
495         &        'BB: not converge: Tsf, dTsf=', Tsf(i,j), dTsrf(i,j)
496               WRITE(stdUnit,*)
497         &        'BB: not converge: flxNet,dFlxT=', flxNet, dFlxdT(i,j)
498               IF ( Tsf(i,j).LT.-70. _d 0 ) THEN
499                 WRITE( msgBuf, '(A,2I4,2I3,I10,F12.3)')
500         &        'THSICE_SOLVE4TEMP: Too low Tsf in', i, j, bi, bj,
501         &                            myIter, Tsf(i,j)
502                 CALL PRINT_ERROR( msgBuf , myThid )
503                 STOP 'ABNORMAL END: S/R THSICE_SOLVE4TEMP'
504               ENDIF
505              ENDIF
506  #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  
507    
508             ENDIF
509            ENDDO
510           ENDDO
511    #ifndef ALLOW_AUTODIFF_TAMC
512         ENDIF         ENDIF
513    #endif
514        ENDDO        ENDDO
515  C ------ end iteration ------------  C ------ end iteration ------------
516    
517  C Compute new bottom layer temperature.  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|
518    
519  #ifdef ALLOW_AUTODIFF_TAMC        DO j = jMin, jMax
520  CADJ STORE  Tice(:)      = comlev1_thsice_1, key=ikey_1         DO i = iMin, iMax
521  CADJ STORE  df0dt        = comlev1_thsice_1, key=ikey_1          IF ( iceMask(i,j).GT.0. _d 0 ) THEN
522  #endif  
523        Tice(2) = (2. _d 0*dt*k32*(Tice(1)+2. _d 0*Tf)  C--   Compute new bottom layer temperature.
524       &        + rhoi*cpIce *hi*Tice(2))            k32       = 2. _d 0*kIce  / hIce(i,j)
525       &         /(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))
526         &                 + rhoi*cpIce*hIce(i,j)*tIc2(i,j))
527         &               /(6. _d 0*dt*k32 + rhoi*cpIce*hIce(i,j))
528  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
529        IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)            IF ( dBug(i,j,bi,bj) ) WRITE(stdUnit,1010)
530       &  'ThSI_SOLVE4T: k, Ts, Tice=',k,Tsf,Tice       &  'ThSI_SOLVE4T: k, Ts, Tice=',k,Tsf(i,j),tIc1(i,j),tIc2(i,j)
531              netSW   = flxAtm(i,j)
532  #endif  #endif
533    
534  C Compute final flux values at surfaces.  C--   Compute final flux values at surfaces.
535              tSrf(i,j)  = Tsf(i,j)
536        fct    = k12*(Tsf-Tice(1))            fct    = k12(i,j)*(Tsf(i,j)-tIc1(i,j))
537        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)
538        flx0   = flx0 + df0dT*dTsf            flxNet = sHeat(i,j) + flxTexSW(i,j)
539        IF ( useBlkFlx ) THEN            flxNet = flxNet + dFlxdT(i,j)*dTsrf(i,j)
540  C--   needs to update also Evap (Tsf changes) since Latent heat has been updated            IF ( useBlkFlx ) THEN
541          evpAtm(i,j) = evap + dEvdT*dTsf  C-    needs to update also Evap (Tsf changes) since Latent heat has been updated
542        ELSE              evpAtm(i,j) = evapT(i,j) + dEvdT(i,j)*dTsrf(i,j)
543              ELSE
544  C- WARNING: Evap & +Evap*Lfresh are missing ! (but only affects Diagnostics)  C- WARNING: Evap & +Evap*Lfresh are missing ! (but only affects Diagnostics)
545          evpAtm(i,j) = 0.              evpAtm(i,j) = 0.
546        ENDIF            ENDIF
547  C-    energy flux to Atmos: use net short-wave flux at surf. and  C-    Update energy flux to Atmos with other than SW contributions;
548  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)
549        flxAtm(i,j) = netSW + flxExceptSw            flxAtm(i,j) = flxAtm(i,j) + flxTexSW(i,j)
550       &                    + df0dT*dTsf + evpAtm(i,j)*Lfresh       &                + dFlxdT(i,j)*dTsrf(i,j) + evpAtm(i,j)*Lfresh
551  C-    excess of energy @ surface (used for surface melting):  C-    excess of energy @ surface (used for surface melting):
552        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  
553    
554  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
555        IF ( dBug(i,j,bi,bj) ) WRITE(6,1020)            IF ( dBug(i,j,bi,bj) ) WRITE(stdUnit,1020)
556       &  'ThSI_SOLVE4T: flx0,fct,Dif,flxCnB=',       &     'ThSI_SOLVE4T: flxNet,fct,Dif,flxCnB=',
557       &                 flx0,fct,flx0-fct,flxCnB(i,j)       &                    flxNet,fct,flxNet-fct,flxCnB(i,j)
558  #endif  #endif
559    
560  C Compute new enthalpy for each layer.  C--   Compute new enthalpy for each layer.
561              qIc1(i,j) = -cpWater*Tmlt1 + cpIce *(Tmlt1-tIc1(i,j))
562         &                + Lfresh*(1. _d 0-Tmlt1/tIc1(i,j))
563              qIc2(i,j) = -cpIce *tIc2(i,j) + Lfresh
564    
       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 i0swFrac)  
   
       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  
565  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
566    C--   Make sure internal ice temperatures do not exceed Tmlt.
567    C     (This should not happen for reasonable values of i0swFrac)
568              IF (tIc1(i,j) .GE. Tmlt1) THEN
569               WRITE(stdUnit,'(A,2I4,2I3,1P2E14.6)')
570         &     ' BBerr - Bug: IceT(1) > Tmlt',i,j,bi,bj,tIc1(i,j),Tmlt1
571              ENDIF
572              IF (tIc2(i,j) .GE. 0. _d 0) THEN
573               WRITE(stdUnit,'(A,2I4,2I3,1P2E14.6)')
574         &     ' BBerr - Bug: IceT(2) > 0',i,j,bi,bj,tIc2(i,j)
575              ENDIF
576    
577            IF ( dBug(i,j,bi,bj) ) THEN            IF ( dBug(i,j,bi,bj) ) THEN
578             WRITE(6,1020) 'ThSI_SOLV_4T: Tsf, Tice(1,2), dTsurf=',             WRITE(stdUnit,1020) 'ThSI_SOLV_4T: Tsf, Tice(1,2), dTsurf=',
579       &                              Tsf, Tice, dTsf       &           Tsf(i,j), tIc1(i,j), tIc2(i,j), dTsrf(i,j)
580             WRITE(6,1020) 'ThSI_SOLV_4T: sHeat, flxCndBt, Qice =',             WRITE(stdUnit,1020) 'ThSI_SOLV_4T: sHeat, flxCndBt, Qice =',
581       &           sHeat(i,j), flxCnB(i,j), qicen       &           sHeat(i,j), flxCnB(i,j), qIc1(i,j), qIc2(i,j)
582             WRITE(6,1020) 'ThSI_SOLV_4T: flxA, evpA, fxSW_bf,af=',             WRITE(stdUnit,1020) 'ThSI_SOLV_4T: flxA, evpA, fxSW_bf,af=',
583       &           flxAtm(i,j), evpAtm(i,j), netSW, flxSW(i,j)       &           flxAtm(i,j), evpAtm(i,j), netSW, flxSW(i,j)
584            ENDIF            ENDIF
585  #endif  #endif
586    
587          ELSE          ELSE
588            IF ( .NOT.useBlkFlx ) dTsrf(i,j) = 0. _d 0  C--     ice-free grid point:
589    c         tIc1  (i,j) = 0. _d 0
590    c         tIc2  (i,j) = 0. _d 0
591              dTsrf (i,j) = 0. _d 0
592    c         sHeat (i,j) = 0. _d 0
593    c         flxCnB(i,j) = 0. _d 0
594    c         flxAtm(i,j) = 0. _d 0
595    c         evpAtm(i,j) = 0. _d 0
596    
597          ENDIF          ENDIF
598         ENDDO         ENDDO
599        ENDDO        ENDDO
600  #endif  /* ALLOW_THSICE */  #endif  /* ALLOW_THSICE */
601    
602  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|
603    
604        RETURN        RETURN
605        END        END
Legend:
Removed from v.1.17  
changed lines
  Added in v.1.27
  ViewVC Help
Powered by ViewVC 1.1.22