/[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.20 by jmc, Thu Sep 30 16:05:19 2010 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,
# 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    C     useBlkFlx    :: use some bulk-formulae to compute fluxes over sea-ice
120    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    C     iceFlag      :: True= do iterate for Surf.Temp ; False= do nothing
123  C     frsnow       :: fractional snow cover  C     frsnow       :: fractional snow cover
124  C     fswpen       :: SW penetrating beneath surface (W m-2)  C     fswpen       :: SW penetrating beneath surface (W m-2)
125  C     fswdn        :: SW absorbed at surface (W m-2)  C     fswdn        :: SW absorbed at surface (W m-2)
126  C     fswint       :: SW absorbed in ice (W m-2)  C     fswint       :: SW absorbed in ice (W m-2)
127  C     fswocn       :: SW passed through ice to ocean (W m-2)  C     fswocn       :: SW passed through ice to ocean (W m-2)
128    C     Tsf          :: surface (ice or snow) temperature   (local copy of tSrf)
129  C     flx0exSW     :: net surface heat flux over melting snow/ice, except short-wave.  C     flx0exSW     :: net surface heat flux over melting snow/ice, except short-wave.
130  C     flxTexSW     :: net surface heat flux, except short-wave (W/m2)  C     flxTexSW     :: net surface heat flux, except short-wave (W/m2)
131    C     dFlxdT       :: deriv of flxNet wrt Tsf (W m-2 deg-1)
132  C     evap0        :: evaporation over melting snow/ice [kg/m2/s] (>0 if evaporate)  C     evap0        :: evaporation over melting snow/ice [kg/m2/s] (>0 if evaporate)
133  C     evapT        :: evaporation over snow/ice [kg/m2/s] (>0 if evaporate)  C     evapT        :: evaporation over snow/ice [kg/m2/s] (>0 if evaporate)
134  C     dEvdT        :: derivative of evap. with respect to Tsf [kg/m2/s/K]  C     dEvdT        :: derivative of evap. with respect to Tsf [kg/m2/s/K]
135  C     flxNet       :: net surf heat flux (+=down), from Atmos. to sea-ice (W m-2)  C     flxNet       :: net surf heat flux (+=down), from Atmos. to sea-ice (W m-2)
136    C     netSW        :: net Short-Wave flux at surface (+=down) [W/m2]
137  C     fct          :: heat conducted to top surface  C     fct          :: heat conducted to top surface
 C     dFlxdT       :: deriv of flxNet wrt Tsf (W m-2 deg-1)  
138  C     k12, k32     :: thermal conductivity terms  C     k12, k32     :: thermal conductivity terms
139  C     a10, b10     :: coefficients in quadratic eqn for T1  C     a10,b10,c10  :: coefficients in quadratic eqn for T1
140  C     a1, b1, c1   :: coefficients in quadratic eqn for T1  C     a1, b1, c1   :: coefficients in quadratic eqn for T1
 C     Tsf_start    :: old value of Tsf  
141  C     dt           :: timestep  C     dt           :: timestep
142          LOGICAL useBlkFlx
143          LOGICAL iterate4Tsf
144          LOGICAL iceFlag(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
145        INTEGER i,j        INTEGER i,j
146        INTEGER k, iterMax        INTEGER k, iterMax, ii, jj, icount
147        _RL  frsnow        _RL  frsnow
148        _RL  fswpen        _RL  fswpen
149        _RL  fswdn        _RL  fswdn
150        _RL  fswint        _RL  fswint
151        _RL  fswocn        _RL  fswocn
152        _RL  flx0exSW        _RL  Tsf     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
153        _RL  flxTexSW        _RL  flx0exSW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
154        _RL  evap0, evapT, dEvdT        _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        _RL  flxNet
160        _RL  fct        _RL  fct
161        _RL  dFlxdT        _RL  k12     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
162        _RL  k12, k32        _RL  k32
163        _RL  a10, b10        _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        _RL  a1, b1, c1
 c     _RL  Tsf_start  
167        _RL  dt        _RL  dt
168        _RL  recip_dhSnowLin        _RL  recip_dhSnowLin
169        LOGICAL useBlkFlx  #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    #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        useBlkFlx = useEXF .OR. useBulkForce
205          dt  = thSIce_dtTemp
206        IF ( dhSnowLin.GT.0. _d 0 ) THEN        IF ( dhSnowLin.GT.0. _d 0 ) THEN
207          recip_dhSnowLin = 1. _d 0 / dhSnowLin          recip_dhSnowLin = 1. _d 0 / dhSnowLin
208        ELSE        ELSE
209          recip_dhSnowLin = 0. _d 0          recip_dhSnowLin = 0. _d 0
210        ENDIF        ENDIF
211    
212        dt  = thSIce_dtTemp        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  dFlxdT        = 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)  
           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.hIceMin ) THEN            IF ( hIce(i,j).LT.hIceMin ) THEN
242  C If hi < hIceMin, melt the ice.  C     if hi < hIceMin, melt the ice.
243               STOP 'THSICE_SOLVE4TEMP: should not enter if hi<hIceMin'  C     keep the position of this problem but do the stop later
244               ii = i
245               jj = j
246               icount = icount + 1
247            ENDIF            ENDIF
248    
249  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  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
259    
260  C fractional snow cover  C--   Compute SW flux absorbed at surface and penetrating to layer 1.
261  C  assume a linear distribution of snow thickness, with dhSnowLin slope,            fswpen = flxSW(i,j) * (1. _d 0 - frsnow) * i0swFrac
262  C   from hs-dhSnowLin to hs+dhSnowLin if full ice & snow cover.            fswocn = fswpen * exp(-ksolar*hIce(i,j))
263  C frsnow = fraction of snow over the ice-covered part of the grid cell            fswint = fswpen - fswocn
264        IF ( hs .GT. iceMask(i,j)*dhSnowLin ) THEN            fswdn  = flxSW(i,j) - fswpen
265          frsnow = 1. _d 0  
266        ELSE  C     Initialise Atmospheric surf. heat flux with net SW flux at the surface
267          frsnow = hs*recip_dhSnowLin/iceMask(i,j)            flxAtm(i,j) = flxSW(i,j)
268          IF ( frsnow.GT.0. _d 0 ) frsnow = SQRT(frsnow)  C     SW flux at sea-ice base left to the ocean
269        ENDIF            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 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  
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,tSrf(i,j),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              c10(i,j) = rhoi*Lfresh*hIce(i,j)*Tmlt1 / (2. _d 0*dt)
315    
316  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|          ELSE
317  C get surface fluxes over melting surface            iceFlag(i,j) = .FALSE.
318            IF ( useBlkFlx ) THEN          ENDIF
319             Tsf = 0.         ENDDO
320             IF ( useEXF ) THEN        ENDDO
321              CALL THSICE_GET_EXF(        IF ( icount .gt. 0 ) THEN
322       I                            hSnow(i,j), Tsf,         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
330    
331    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|
332    
333    #ifdef ALLOW_AUTODIFF_TAMC
334    CADJ STORE devdt    = comlev1_bibj,key=iicekey,byte=isbyte
335    CADJ STORE tsf      = comlev1_bibj,key=iicekey,byte=isbyte
336    #endif
337    
338    C--   Get surface fluxes over melting surface
339    #ifdef ALLOW_AUTODIFF_TAMC
340          IF ( useBlkFlx ) THEN
341    #else
342          IF ( useBlkFlx .AND. iterate4Tsf  ) THEN
343    #endif
344            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,       O                            flx0exSW, dFlxdT, evap0, dEvdT,
354       I                            i,j,bi,bj,myThid )       I                            myTime, myIter, myThid )
355  C could add this "ifdef" to hide THSICE_GET_BULKF from TAF  C-    could add this "ifdef" to hide THSICE_GET_BULKF from TAF
356  c#ifdef ALLOW_BULK_FORCE  #ifdef ALLOW_BULK_FORCE
357             ELSEIF ( useBulkForce ) THEN          ELSEIF ( useBulkForce ) THEN
358              CALL THSICE_GET_BULKF(             CALL THSICE_GET_BULKF( bi, bj,
359       I                            hSnow(i,j), Tsf,       I                            iMin, iMax, jMin, jMax,
360         I                            iceFlag, hSnow, Tsf,
361       O                            flx0exSW, dFlxdT, evap0, dEvdT,       O                            flx0exSW, dFlxdT, evap0, dEvdT,
362       I                            i,j,bi,bj,myThid )       I                            myTime, myIter, myThid )
363  c#endif /* ALLOW_BULK_FORCE */  #endif /* ALLOW_BULK_FORCE */
364             ENDIF          ENDIF
365            ELSE        ENDIF
            flx0exSW = flxExSW(i,j,0)  
           ENDIF  
   
 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  
 C Compute new surface and internal temperatures; iterate until  
 C Tsfc converges.  
366    
367    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|
368    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        IF ( useBlkFlx ) THEN
377          iterMax = nitMaxTsf          iterMax = nitMaxTsf
378        ELSE        ELSE
379          iterMax = 1          iterMax = 1
380        ENDIF        ENDIF
       Tsf  = tSrf(i,j)  
       dTsf = Terrmax  
   
 C ----- begin iteration  -----  
       DO k = 1,iterMax  
381    
382  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
383         ikey_3 = (ikey_1-1)*MaxTsf + k  CADJ STORE devdt    = comlev1_bibj,key=iicekey,byte=isbyte
384    CADJ STORE dflxdt   = comlev1_bibj,key=iicekey,byte=isbyte
385    CADJ STORE flx0exsw = comlev1_bibj,key=iicekey,byte=isbyte
386    CADJ STORE flxtexsw = comlev1_bibj,key=iicekey,byte=isbyte
387  #endif  #endif
388    
389  #ifdef ALLOW_AUTODIFF_TAMC  C ----- begin iteration  -----
390  CADJ STORE tsf          = comlev1_thsice_3, key=ikey_3        DO k = 1,iterMax
391  CADJ STORE dtsf         = comlev1_thsice_3, key=ikey_3  #ifndef ALLOW_AUTODIFF_TAMC
392  CADJ STORE dFlxdT        = comlev1_thsice_3, key=ikey_3         IF ( iterate4Tsf ) THEN
393  CADJ STORE flxexceptsw  = comlev1_thsice_3, key=ikey_3         iterate4Tsf = .FALSE.
394  #endif  #endif
        IF ( ABS(dTsf).GE.Terrmax ) THEN  
395    
396  C Save temperatures at start of iteration.         IF ( useBlkFlx ) THEN
397  c        Tsf_start = Tsf  C--   Compute top surface flux.
   
         IF ( useBlkFlx ) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE tsf          = comlev1_thsice_3, key=ikey_3  
 #endif  
 C Compute top surface flux.  
398           IF ( useEXF ) THEN           IF ( useEXF ) THEN
399            CALL THSICE_GET_EXF  (             CALL THSICE_GET_EXF(   bi, bj,
400       I                         hs, Tsf,       I                            iMin, iMax, jMin, jMax,
401       O                         flxTexSW, dFlxdT, evapT, dEvdT,       I                            iceFlag, hSnow, Tsf,
402       I                         i,j,bi,bj,myThid )       O                            flxTexSW, dFlxdT, evapT, dEvdT,
403  C could add this "ifdef" to hide THSICE_GET_BULKF from TAF       I                            myTime, myIter, myThid )
404  c#ifdef ALLOW_BULK_FORCE  C-    could add this "ifdef" to hide THSICE_GET_BULKF from TAF
405    #ifdef ALLOW_BULK_FORCE
406           ELSEIF ( useBulkForce ) THEN           ELSEIF ( useBulkForce ) THEN
407            CALL THSICE_GET_BULKF(             CALL THSICE_GET_BULKF( bi, bj,
408       I                         hs, Tsf,       I                            iMin, iMax, jMin, jMax,
409       O                         flxTexSW, dFlxdT, evapT, dEvdT,       I                            iceFlag, hSnow, Tsf,
410       I                         i,j,bi,bj,myThid )       O                            flxTexSW, dFlxdT, evapT, dEvdT,
411  c#endif /* ALLOW_BULK_FORCE */       I                            myTime, myIter, myThid )
412    #endif /* ALLOW_BULK_FORCE */
413           ENDIF           ENDIF
414          ELSE         ELSE
415           flxTexSW = flxExSW(i,j,1)           DO j = jMin, jMax
416           dFlxdT = flxExSW(i,j,2)            DO i = iMin, iMax
417          ENDIF             IF ( iceFlag(i,j) ) THEN
418          flxNet = fswdn + flxTexSW               flxTexSW(i,j) = flxExSW(i,j,1)
419  #ifdef ALLOW_DBUG_THSICE               dFlxdT(i,j) = flxExSW(i,j,2)
420        IF ( dBug(i,j,bi,bj) ) WRITE(6,1020)             ENDIF
421       &  'ThSI_SOLVE4T: flxNet,dFlxdT,k12,D=',            ENDDO
422       &                 flxNet, dFlxdT, k12, k12-dFlxdT           ENDDO
423  #endif         ENDIF
   
 C Compute new top layer and surface temperatures.  
 C If Tsfc is computed to be > 0 C, fix Tsfc = 0 and recompute T1  
 C with different coefficients.  
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
432    
433    C--   Compute new top layer and surface temperatures.
434    C     If Tsfc is computed to be > 0 C, fix Tsfc = 0 and recompute T1
435    C     with different coefficients.
436           DO j = jMin, jMax
437            DO i = iMin, iMax
438             IF ( iceFlag(i,j) ) THEN
439              flxNet = sHeat(i,j) + flxTexSW(i,j)
440    #ifdef ALLOW_DBUG_THSICE
441              IF ( dBug(i,j,bi,bj) ) WRITE(6,1020)
442         &     '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*dFlxdT / (k12-dFlxdT)  
446           b1 = b10 - k12*(flxNet-dFlxdT*Tsf) / (k12-dFlxdT)            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           dTsf = (flxNet + k12*(Tice(1)-Tsf)) / (k12-dFlxdT)       &                  /(k12(i,j)-dFlxdT(i,j))
449           Tsf = Tsf + dTsf            c1 = c10(i,j)
450           IF (Tsf .GT. 0. _d 0) THEN            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  #ifdef ALLOW_DBUG_THSICE
456              IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)             IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)
457       &        '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)
458  #endif  #endif
459              a1 = a10 + k12             a1 = a10(i,j) + k12(i,j)
460  C      note: b1 = b10 - k12*Tf0  C      note: b1 = b10 - k12*Tf0
461              b1 = b10             b1 = b10(i,j)
462              Tice(1) = (-b1 - SQRT(b1*b1-4. _d 0*a1*c1))/(2. _d 0*a1)             tIc1(i,j) = (-b1 - SQRT(b1*b1-4. _d 0*a1*c1))/(2. _d 0*a1)
463              Tsf = 0. _d 0             Tsf(i,j) = 0. _d 0
464             IF ( useBlkFlx ) THEN             IF ( useBlkFlx ) THEN
465              flxTexSW = flx0exSW              flxTexSW(i,j) = flx0exSW(i,j)
466              evapT = evap0              evapT(i,j) = evap0(i,j)
467              dTsf = 0. _d 0              dTsrf(i,j) = 0. _d 0
468             ELSE             ELSE
469              flxTexSW = flxExSW(i,j,0)              flxTexSW(i,j) = flxExSW(i,j,0)
470              dTsf = 1000.              dTsrf(i,j) = 1000.
471              dFlxdT = 0.              dFlxdT(i,j) = 0.
472             ENDIF             ENDIF
473             flxNet = fswdn + flxTexSW            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: flxNet,dFlxT=',flxNet,dFlxdT  
             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  dFlxdT        = 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        flxNet   = flxNet + dFlxdT*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) = evapT + 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 + flxTexSW            flxAtm(i,j) = flxAtm(i,j) + flxTexSW(i,j)
540       &                    + dFlxdT*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) = flxNet - 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: flxNet,fct,Dif,flxCnB=',       &     'ThSI_SOLVE4T: flxNet,fct,Dif,flxCnB=',
547       &                 flxNet,fct,flxNet-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        qicen(1) = -cpWater*Tmlt1 + cpIce *(Tmlt1-Tice(1))       &                + Lfresh*(1. _d 0-Tmlt1/tIc1(i,j))
553       &            + Lfresh*(1. _d 0-Tmlt1/Tice(1))            qIc2(i,j) = -cpIce *tIc2(i,j) + Lfresh
       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  
554    
 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.20  
changed lines
  Added in v.1.24
  ViewVC Help
Powered by ViewVC 1.1.22