/[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.7 by jmc, Mon Mar 13 03:55:39 2006 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                     useBlkFlx, flxExcSw, Tf, hi, hs,       I                  bi, bj,
11       U                     flxSW, Tsf, qicen,       I                  iMin,iMax, jMin,jMax, dBugFlag,
12       O                     Tice, sHeating, flxCnB,       I                  useBulkForce, useEXF,
13       O                     dTsf, flxAtm, evpAtm,       I                  iceMask, hIce, hSnow, tFrz, flxExSW,
14       I                     i,j,bi,bj, myThid)       U                  flxSW, tSrf, qIc1, qIc2,
15         O                  tIc1, tIc2, dTsrf,
16         O                  sHeat, flxCnB, flxAtm, evpAtm,
17         I                  myTime, myIter, myThid )
18  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
19  C     *==========================================================*  C     *==========================================================*
20  C     | S/R  THSICE_SOLVE4TEMP  C     | S/R  THSICE_SOLVE4TEMP
# Line 20  C     | Solve (implicitly) for sea-ice a Line 23  C     | Solve (implicitly) for sea-ice a
23  C     *==========================================================*  C     *==========================================================*
24  C     \ev  C     \ev
25    
 C     !USES:  
       IMPLICIT NONE  
   
 C     == Global variables ===  
 #include "EEPARAMS.h"  
 #include "THSICE_SIZE.h"  
 #include "THSICE_PARAMS.h"  
   
 C     !INPUT/OUTPUT PARAMETERS:  
 C     == Routine Arguments ==  
 C    useBlkFlx :: use surf. fluxes from bulk-forcing external S/R  
 C     flxExcSw :: surf. heat flux (+=down) except SW, function of surf. temp Ts:  
 C                0: Flx(Ts=0) ; 1: Flx(Ts=Ts^n) ; 2: d.Flx/dTs(Ts=Ts^n)  
 C     Tf       :: freezing temperature (oC) of local sea-water  
 C     hi       :: ice height  
 C     hs       :: snow height  
 C     flxSW    :: net Short-Wave flux (+=down) [W/m2]: input= at surface  
 C              ::               output= at the sea-ice base to the ocean  
 C     Tsf      :: surface (ice or snow) temperature  
 C     qicen    :: ice enthalpy (J/kg)  
 C     Tice     :: internal ice temperatures  
 C     sHeating :: surf heating left to melt snow or ice (= Atmos-conduction)  
 C     flxCnB   :: heat flux conducted through the ice to bottom surface  
 C     dTsf     :: surf. temp adjusment: Ts^n+1 - Ts^n  
 C     flxAtm   :: net flux of energy from the atmosphere [W/m2] (+=down)  
 C                without snow precip. (energy=0 for liquid water at 0.oC)  
 C     evpAtm   :: evaporation to the atmosphere (kg/m2/s) (>0 if evaporate)  
 C   i,j,bi,bj  :: indices of current grid point  
 C     myThid   :: Thread no. that called this routine.  
       LOGICAL useBlkFlx  
       _RL flxExcSw(0:2)  
       _RL Tf  
       _RL hi  
       _RL hs  
   
       _RL flxSW  
       _RL Tsf  
       _RL qicen(nlyr)  
   
       _RL Tice (nlyr)  
       _RL sHeating  
       _RL flxCnB  
       _RL dTsf  
       _RL flxAtm  
       _RL evpAtm  
       INTEGER i,j, bi,bj  
       INTEGER myThid  
 CEOP  
   
 #ifdef ALLOW_THSICE  
   
26  C ADAPTED FROM:  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 87  Cc****subroutine thermo_winton(n,fice,fs Line 39  Cc****subroutine thermo_winton(n,fice,fs
39  C.. Compute temperature change using Winton model with 2 ice layers, of  C.. Compute temperature change using Winton model with 2 ice layers, of
40  C.. which only the top layer has a variable heat capacity.  C.. which only the top layer has a variable heat capacity.
41    
42  C     == Local Variables ==  C     !USES:
43        INTEGER k, iterMax        IMPLICIT NONE
   
       _RL  frsnow        ! fractional snow cover  
   
       _RL  fswpen        ! SW penetrating beneath surface (W m-2)  
       _RL  fswdn         ! SW absorbed at surface (W m-2)  
       _RL  fswint        ! SW absorbed in ice (W m-2)  
       _RL  fswocn        ! SW passed through ice to ocean (W m-2)  
   
       _RL  flxExceptSw   ! net surface heat flux, except short-wave (W/m2)  
 C     evap           :: evaporation over snow/ice [kg/m2/s] (>0 if evaporate)  
 C     dEvdT          :: derivative of evap. with respect to Tsf [kg/m2/s/K]  
       _RL  evap, dEvdT  
       _RL  flx0          ! net surf heat flux, from Atmos. to sea-ice (W m-2)  
       _RL  fct           ! heat conducted to top surface  
44    
45        _RL  df0dT         ! deriv of flx0 wrt Tsf (W m-2 deg-1)  C     == Global variables ===
46    #include "EEPARAMS.h"
47    #include "SIZE.h"
48    #include "THSICE_SIZE.h"
49    #include "THSICE_PARAMS.h"
50    #ifdef ALLOW_AUTODIFF_TAMC
51    # include "tamc.h"
52    # include "tamc_keys.h"
53    #endif
54    
55        _RL  k12, k32      ! thermal conductivity terms  C     !INPUT/OUTPUT PARAMETERS:
56        _RL  a10, b10      ! coefficients in quadratic eqn for T1  C     == Routine Arguments ==
57        _RL  a1, b1, c1    ! coefficients in quadratic eqn for T1  C     bi,bj       :: tile indices
58  c     _RL  Tsf_start     ! old value of Tsf  C     iMin,iMax   :: computation domain: 1rst index range
59    C     jMin,jMax   :: computation domain: 2nd  index range
60    C     dBugFlag    :: allow to print debugging stuff (e.g. on 1 grid point).
61    C     useBulkForce:: use surf. fluxes from bulk-forcing external S/R
62    C     useEXF      :: use surf. fluxes from exf          external S/R
63    C---  Input:
64    C     iceMask     :: sea-ice fractional mask [0-1]
65    C     hIce        :: ice height [m]
66    C     hSnow       :: snow height [m]
67    C     tFrz        :: sea-water freezing temperature [oC] (function of S)
68    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)
70    C---  Modified (input&output):
71    C     flxSW       :: net Short-Wave flux (+=down) [W/m2]: input= at surface
72    C                 ::                 output= below sea-ice, into the ocean
73    C     tSrf        :: surface (ice or snow) temperature
74    C     qIc1        :: ice enthalpy (J/kg), 1rst level
75    C     qIc2        :: ice enthalpy (J/kg), 2nd level
76    C---  Output
77    C     tIc1        :: temperature of ice layer 1 [oC]
78    C     tIc2        :: temperature of ice layer 2 [oC]
79    C     dTsrf       :: surf. temp adjusment: Ts^n+1 - Ts^n
80    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
82    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)
84    C     evpAtm      :: evaporation to the atmosphere [kg/m2/s] (>0 if evaporate)
85    C---  Input:
86    C     myTime      :: current Time of simulation [s]
87    C     myIter      :: current Iteration number in simulation
88    C     myThid      :: my Thread Id number
89          INTEGER bi,bj
90          INTEGER iMin, iMax
91          INTEGER jMin, jMax
92          LOGICAL dBugFlag
93          LOGICAL useBulkForce
94          LOGICAL useEXF
95          _RL iceMask(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
96          _RL hIce   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
97          _RL hSnow  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98          _RL tFrz   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
99          _RL flxExSW(iMin:iMax,jMin:jMax,0:2)
100          _RL flxSW  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
101          _RL tSrf   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102          _RL qIc1   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103          _RL qIc2   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
104          _RL tIc1   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105          _RL tIc2   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
106          _RL dTsrf  (iMin:iMax,jMin:jMax)
107          _RL sHeat  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
108          _RL flxCnB (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
109          _RL flxAtm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
110          _RL evpAtm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
111          _RL  myTime
112          INTEGER myIter
113          INTEGER myThid
114    CEOP
115    
116        _RL  dt            ! timestep  #ifdef ALLOW_THSICE
117    C     !LOCAL VARIABLES:
118    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
124    C     fswpen       :: SW penetrating beneath surface (W m-2)
125    C     fswdn        :: SW absorbed at surface (W m-2)
126    C     fswint       :: SW absorbed in ice (W m-2)
127    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.
130    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)
133    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]
135    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
138    C     k12, k32     :: thermal conductivity terms
139    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
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        INTEGER iceornot  C-    Define grid-point location where to print debugging values
174        LOGICAL dBug  #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        dt  = thSIce_deltaT  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|
180        dBug = .FALSE.  
181  c     dBug = ( bi.EQ.3 .AND. i.EQ.15 .AND. j.EQ.11 )  #ifdef ALLOW_AUTODIFF_TAMC
182  c     dBug = ( bi.EQ.6 .AND. i.EQ.10 .AND. j.EQ.20 )  c     act1 = bi - myBxLo(myThid)
183        IF (dBug) WRITE(6,'(A,2I4,2I2)') 'ThSI_SOLVE4T: i,j=',i,j,bi,bj  c     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
184    c     act2 = bj - myByLo(myThid)
185        IF ( hi.LT.himin ) THEN  c     max2 = myByHi(myThid) - myByLo(myThid) + 1
186  C If hi < himin, melt the ice.  c     act3 = myThid - 1
187           STOP 'THSICE_SOLVE4TEMP: should not enter if hi<himin'  c     max3 = nTx*nTy
188    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 */
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
205          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        ENDIF
211    
212  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|        iterate4Tsf = .FALSE.
213          icount = 0
214    C    
215          DO j = jMin, jMax
216           DO i = iMin, iMax
217    #ifdef ALLOW_AUTODIFF_TAMC
218    c         ikey_1 = i
219    c    &         + sNx*(j-1)
220    c    &         + sNx*sNy*act1
221    c    &         + sNx*sNy*max1*act2
222    c    &         + sNx*sNy*max1*max2*act3
223    c    &         + sNx*sNy*max1*max2*max3*act4
224    #endif /* ALLOW_AUTODIFF_TAMC */
225    #ifdef ALLOW_AUTODIFF_TAMC
226    cCADJ STORE  devdt        = comlev1_thsice_1, key=ikey_1
227    cCADJ STORE  dFlxdT        = comlev1_thsice_1, key=ikey_1
228    cCADJ STORE  flxexceptsw  = comlev1_thsice_1, key=ikey_1
229    cCADJ STORE  flxsw(i,j)   = comlev1_thsice_1, key=ikey_1
230    cCADJ STORE  qic1(i,j)    = comlev1_thsice_1, key=ikey_1
231    cCADJ STORE  qic2(i,j)    = comlev1_thsice_1, key=ikey_1
232    cCADJ STORE  tsrf(i,j)    = comlev1_thsice_1, key=ikey_1
233    #endif
234            IF ( iceMask(i,j).GT.0. _d 0) THEN
235              iterate4Tsf  = .TRUE.
236              iceFlag(i,j) = .TRUE.
237    #ifdef ALLOW_DBUG_THSICE
238              IF ( dBug(i,j,bi,bj) ) WRITE(6,'(A,2I4,2I2)')
239         &         'ThSI_SOLVE4T: i,j=',i,j,bi,bj
240    #endif
241              IF ( hIce(i,j).LT.hIceMin ) THEN
242    C     if hi < hIceMin, melt the ice.
243    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
259    
260    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    
285    #ifdef ALLOW_DBUG_THSICE
286              IF (tIc1(i,j).GT.0. _d 0 ) THEN
287               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
301    
302    C--   Compute coefficients used in quadratic formula.
303    
304              a10(i,j) = rhoi*cpIce *hIce(i,j)/(2. _d 0*dt) +
305         &          k32*( 4. _d 0*dt*k32 + rhoi*cpIce *hIce(i,j) )
306         &           / ( 6. _d 0*dt*k32 + rhoi*cpIce *hIce(i,j) )
307              b10(i,j) = -hIce(i,j)*
308         &          ( rhoi*cpIce*tIc1(i,j) + rhoi*Lfresh*Tmlt1/tIc1(i,j) )
309         &           /(2. _d 0*dt)
310         &        - k32*( 4. _d 0*dt*k32*tFrz(i,j)
311         &               +rhoi*cpIce*hIce(i,j)*tIc2(i,j) )
312         &           / ( 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 fractional snow cover          ELSE
317        frsnow = 0. _d 0            iceFlag(i,j) = .FALSE.
318        IF (hs .GT. 0. _d 0) frsnow = 1. _d 0          ENDIF
319           ENDDO
320  C Compute SW flux absorbed at surface and penetrating to layer 1.        ENDDO
321        fswpen  = flxSW * (1. _d 0 - frsnow) * i0        IF ( icount .gt. 0 ) THEN
322        fswocn = fswpen * exp(-ksolar*hi)         WRITE (standardMessageUnit,'(A,I5,A)')
323        fswint = fswpen - fswocn       &      'THSICE_SOLVE4TEMP: there are ',icount,
324         &      ' case(s) where hIce<hIceMin;'
325        fswdn = flxSW - fswpen         WRITE (standardMessageUnit,'(A,I3,A1,I3,A)')
326         &      'THSICE_SOLVE4TEMP: the last one was at (',ii,',',jj,
327  C Compute conductivity terms at layer interfaces.       &      ') with hIce = ', hIce(ii,jj)
328           STOP 'THSICE_SOLVE4TEMP: should not enter if hIce<hIceMin'
       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 .OR. Tice(2).GT.0. _d 0) THEN  
           WRITE (6,*) 'BBerr Tice(1) > 0 = ',Tice(1)  
           WRITE (6,*) 'BBerr Tice(2) > 0 = ',Tice(2)  
329        ENDIF        ENDIF
       IF (dBug) WRITE(6,1010) 'ThSI_SOLVE4T: k, Ts, Tice=',0,Tsf,Tice  
330    
331  C Compute coefficients used in quadratic formula.  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|
332    
333        a10 = rhoi*cpice *hi/(2. _d 0*dt) +  #ifdef ALLOW_AUTODIFF_TAMC
334       &      k32 * (4. _d 0*dt*k32 + rhoi*cpice *hi)  CADJ STORE devdt    = comlev1_bibj,key=iicekey,byte=isbyte
335       &       / (6. _d 0*dt*k32 + rhoi*cpice *hi)  CADJ STORE tsf      = comlev1_bibj,key=iicekey,byte=isbyte
336        b10 = -hi*  #endif
      &      (rhoi*cpice*Tice(1)+rhoi*Lfresh*Tmlt1/Tice(1))  
      &       /(2. _d 0*dt)  
      &      - k32 * (4. _d 0*dt*k32*Tf+rhoi*cpice *hi*Tice(2))  
      &       / (6. _d 0*dt*k32 + rhoi*cpice *hi) - fswint  
       c1 = rhoi*Lfresh*hi*Tmlt1 / (2. _d 0*dt)  
   
 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  
 C Compute new surface and internal temperatures; iterate until  
 C Tsfc converges.  
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,
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---+----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
381        dTsf = Terrmax  
382    #ifdef ALLOW_AUTODIFF_TAMC
383    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
388    
389  C ----- begin iteration  -----  C ----- begin iteration  -----
390        DO k = 1,iterMax        DO k = 1,iterMax
391         IF ( ABS(dTsf).GE.Terrmax ) THEN  #ifndef ALLOW_AUTODIFF_TAMC
392           IF ( iterate4Tsf ) THEN
393  C Save temperatures at start of iteration.         iterate4Tsf = .FALSE.
394  c        Tsf_start = Tsf  #endif
395    
396          IF ( useBlkFlx ) THEN         IF ( useBlkFlx ) THEN
397  C Compute top surface flux.  C--   Compute top surface flux.
398           IF (hs.GT.3. _d -1) THEN           IF ( useEXF ) THEN
399                iceornot=2             CALL THSICE_GET_EXF(   bi, bj,
400           ELSE       I                            iMin, iMax, jMin, jMax,
401                iceornot=1       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           ENDIF
414           CALL THSICE_GET_BULKF(         ELSE
415       I                        iceornot, Tsf,           DO j = jMin, jMax
416       O                        flxExceptSw, df0dT, evap, dEvdT,            DO i = iMin, iMax
417       I                        i,j,bi,bj,myThid )             IF ( iceFlag(i,j) ) THEN
418          ELSE               flxTexSW(i,j) = flxExSW(i,j,1)
419           flxExceptSw = flxExcSw(1)               dFlxdT(i,j) = flxExSW(i,j,2)
420           df0dT = flxExcSw(2)             ENDIF
421          ENDIF            ENDDO
422          flx0 = fswdn + flxExceptSw           ENDDO
423          IF ( dBug ) WRITE(6,1020) 'ThSI_SOLVE4T: flx0,df0dT,k12,D=',         ENDIF
424       &                                     flx0,df0dT,k12,k12-df0dT  
425    #ifdef ALLOW_AUTODIFF_TAMC
426  C Compute new top layer and surface temperatures.  CADJ STORE devdt(i,j) = comlev1_bibj,key=iicekey,byte=isbyte
427  C If Tsfc is computed to be > 0 C, fix Tsfc = 0 and recompute T1  CADJ STORE dflxdt(i,j) = comlev1_bibj,key=iicekey,byte=isbyte
428  C with different coefficients.  CADJ STORE flxtexsw(i,j) = comlev1_bibj,key=iicekey,byte=isbyte
429    CADJ STORE iceflag(i,j) = comlev1_bibj,key=iicekey,byte=isbyte
430           a1 = a10 - k12*df0dT / (k12-df0dT)  CADJ STORE tsf(i,j) = comlev1_bibj,key=iicekey,byte=isbyte
431           b1 = b10 - k12*(flx0-df0dT*Tsf) / (k12-df0dT)  #endif
432           Tice(1) = -(b1 + SQRT(b1*b1-4. _d 0*a1*c1))/(2. _d 0*a1)  
433           dTsf = (flx0 + k12*(Tice(1)-Tsf)) / (k12-df0dT)  C--   Compute new top layer and surface temperatures.
434           Tsf = Tsf + dTsf  C     If Tsfc is computed to be > 0 C, fix Tsfc = 0 and recompute T1
435           IF (Tsf .GT. 0. _d 0) THEN  C     with different coefficients.
436              IF(dBug) WRITE(6,1010) 'ThSI_SOLVE4T: k,ts,t1,dTs=',         DO j = jMin, jMax
437       &                                          k,Tsf,Tice(1),dTsf          DO i = iMin, iMax
438              a1 = a10 + k12           IF ( iceFlag(i,j) ) THEN
439              b1 = b10          ! note b1 = b10 - k12*Tf0            flxNet = sHeat(i,j) + flxTexSW(i,j)
440              Tice(1) = (-b1 - SQRT(b1*b1-4. _d 0*a1*c1))/(2. _d 0*a1)  #ifdef ALLOW_DBUG_THSICE
441              Tsf = 0. _d 0            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
445    
446              a1 = a10(i,j) - k12(i,j)*dFlxdT(i,j) / (k12(i,j)-dFlxdT(i,j))
447              b1 = b10(i,j) - k12(i,j)*(flxNet-dFlxdT(i,j)*Tsf(i,j))
448         &                  /(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  
             CALL THSICE_GET_BULKF(  
      I                        iceornot, Tsf,  
      O                        flxExceptSw, df0dT, evap, dEvdT,  
      I                        i,j,bi,bj,myThid )  
             dTsf = 0. _d 0  
468             ELSE             ELSE
469              flxExceptSw = flxExcSw(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           IF ( dBug ) WRITE(6,1010) 'ThSI_SOLVE4T: k,ts,t1,dTs=',  C     If no convergence, then repeat.
483       &                                            k,Tsf,Tice(1),dTsf  
484           IF ( useBlkFlx .AND. k.EQ.nitMaxTsf  #ifdef ALLOW_DBUG_THSICE
485       &                  .AND. ABS(dTsf).GE.Terrmax ) THEN            IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)
486              WRITE (6,*) 'BB: thermw conv err ',i,j,bi,bj,dTsf       &    'ThSI_SOLVE4T: k,ts,t1,dTs=', k,Tsf(i,j),tIc1(i,j),dTsrf(i,j)
487              WRITE (6,*) 'BB: thermw conv err, iceheight ', hi            IF ( useBlkFlx .AND. k.EQ.nitMaxTsf .AND. iceFlag(i,j) ) THEN
488              WRITE (6,*) 'BB: thermw conv err: Tsf, flx0', Tsf,flx0             WRITE (6,'(A,4I4,I12,F15.9)')
489              IF (Tsf.LT.-70. _d 0) STOP       &       ' BB: not converge: i,j,it,hi=',i,j,bi,bj,myIter,hIce(i,j)
490           ENDIF             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
497    
498  100   continue  ! surface temperature iteration           ENDIF
499            ENDDO
500           ENDDO
501    #ifndef ALLOW_AUTODIFF_TAMC
502         ENDIF         ENDIF
503    #endif
504        ENDDO        ENDDO
 150   continue  
505  C ------ end iteration ------------  C ------ end iteration ------------
506    
507  C Compute new bottom layer temperature.  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|
   
       Tice(2) = (2. _d 0*dt*k32*(Tice(1)+2. _d 0*Tf)  
      &        + rhoi*cpice *hi*Tice(2))  
      &         /(6. _d 0*dt*k32 + rhoi*cpice *hi)  
       IF (dBug) WRITE(6,1010) 'ThSI_SOLVE4T: k, Ts, Tice=',k,Tsf,Tice  
   
   
 C Compute final flux values at surfaces.  
508    
509        fct    = k12*(Tsf-Tice(1))        DO j = jMin, jMax
510        flxCnB = 4. _d 0*kice *(Tice(2)-Tf)/hi         DO i = iMin, iMax
511        flx0   = flx0 + df0dT*dTsf          IF ( iceMask(i,j).GT.0. _d 0 ) THEN
512        IF ( useBlkFlx ) THEN  
513  C--   needs to update also Evap (Tsf changes) since Latent heat has been updated  C--   Compute new bottom layer temperature.
514          evpAtm = evap + dEvdT*dTsf            k32       = 2. _d 0*kIce  / hIce(i,j)
515        ELSE            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
519              IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)
520         &  'ThSI_SOLVE4T: k, Ts, Tice=',k,Tsf(i,j),tIc1(i,j),tIc2(i,j)
521              netSW   = flxAtm(i,j)
522    #endif
523    
524    C--   Compute final flux values at surfaces.
525              tSrf(i,j)  = Tsf(i,j)
526              fct    = k12(i,j)*(Tsf(i,j)-tIc1(i,j))
527              flxCnB(i,j) = 4. _d 0*kIce *(tIc2(i,j)-tFrz(i,j))/hIce(i,j)
528              flxNet = sHeat(i,j) + flxTexSW(i,j)
529              flxNet = flxNet + dFlxdT(i,j)*dTsrf(i,j)
530              IF ( useBlkFlx ) THEN
531    C-    needs to update also Evap (Tsf changes) since Latent heat has been updated
532                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 = 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 = flxSW + flxExceptSw + df0dT*dTsf + evpAtm*Lfresh            flxAtm(i,j) = flxAtm(i,j) + flxTexSW(i,j)
540         &                + 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        sHeating = flx0 - fct            sHeat(i,j) = flxNet - fct
   
 C-    SW flux at sea-ice base left to the ocean  
       flxSW = fswocn  
543    
544        IF (dBug) WRITE(6,1020) 'ThSI_SOLVE4T: flx0,fct,Dif,flxCnB=',  #ifdef ALLOW_DBUG_THSICE
545       &    flx0,fct,flx0-fct,flxCnB            IF ( dBug(i,j,bi,bj) ) WRITE(6,1020)
546         &     'ThSI_SOLVE4T: flxNet,fct,Dif,flxCnB=',
547         &                    flxNet,fct,flxNet-fct,flxCnB(i,j)
548    #endif
549    
550    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    
555    #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
568               WRITE(6,1020) 'ThSI_SOLV_4T: Tsf, Tice(1,2), dTsurf=',
569         &           Tsf(i,j), tIc1(i,j), tIc2(i,j), dTsrf(i,j)
570               WRITE(6,1020) 'ThSI_SOLV_4T: sHeat, flxCndBt, Qice =',
571         &           sHeat(i,j), flxCnB(i,j), qIc1(i,j), qIc2(i,j)
572               WRITE(6,1020) 'ThSI_SOLV_4T: flxA, evpA, fxSW_bf,af=',
573         &           flxAtm(i,j), evpAtm(i,j), netSW, flxSW(i,j)
574              ENDIF
575    #endif
576    
577  C Compute new enthalpy for each layer.          ELSE
578    C--     ice-free grid point:
579        qicen(1) = -cpwater*Tmlt1 + cpice *(Tmlt1-Tice(1))  c         tIc1  (i,j) = 0. _d 0
580       &            + Lfresh*(1. _d 0-Tmlt1/Tice(1))  c         tIc2  (i,j) = 0. _d 0
581        qicen(2) = -cpice *Tice(2) + Lfresh            dTsrf (i,j) = 0. _d 0
582    c         sHeat (i,j) = 0. _d 0
583  C Make sure internal ice temperatures do not exceed Tmlt.  c         flxCnB(i,j) = 0. _d 0
584  C (This should not happen for reasonable values of i0.)  c         flxAtm(i,j) = 0. _d 0
585    c         evpAtm(i,j) = 0. _d 0
       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  
586    
587            ENDIF
588           ENDDO
589          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.7  
changed lines
  Added in v.1.24
  ViewVC Help
Powered by ViewVC 1.1.22