/[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.5 by jmc, Fri Dec 17 03:44:52 2004 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    
26    C ADAPTED FROM:
27    C LANL CICE.v2.0.2
28    C-----------------------------------------------------------------------
29    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
31    C..       thermodynamic sea ice model for climate study.
32    C..       J. Geophys. Res., 104, 15669 - 15677.
33    C..     Winton, M., 1999:  "A reformulated three-layer sea ice model."
34    C..       Submitted to J. Atmos. Ocean. Technol.
35    C.. authors Elizabeth C. Hunke and William Lipscomb
36    C..         Fluid Dynamics Group, Los Alamos National Laboratory
37    C-----------------------------------------------------------------------
38    Cc****subroutine thermo_winton(n,fice,fsnow,dqice,dTsfc)
39    C.. Compute temperature change using Winton model with 2 ice layers, of
40    C.. which only the top layer has a variable heat capacity.
41    
42  C     !USES:  C     !USES:
43        IMPLICIT NONE        IMPLICIT NONE
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
51    # include "tamc.h"
52    # include "tamc_keys.h"
53    #endif
54    
55  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
56  C     == Routine Arguments ==  C     == Routine Arguments ==
57  C    useBlkFlx :: use surf. fluxes from bulk-forcing external S/R  C     bi,bj       :: tile indices
58  C     flxExcSw :: surf. heat flux (+=down) except SW, function of surf. temp Ts:  C     iMin,iMax   :: computation domain: 1rst index range
59  C                0: Flx(Ts=0) ; 1: Flx(Ts=Ts^n) ; 2: d.Flx/dTs(Ts=Ts^n)  C     jMin,jMax   :: computation domain: 2nd  index range
60  C     Tf       :: freezing temperature (oC) of local sea-water  C     dBugFlag    :: allow to print debugging stuff (e.g. on 1 grid point).
61  C     hi       :: ice height  C     useBulkForce:: use surf. fluxes from bulk-forcing external S/R
62  C     hs       :: snow height  C     useEXF      :: use surf. fluxes from exf          external S/R
63  C     flxSW    :: net Short-Wave flux (+=down) [W/m2]: input= at surface  C---  Input:
64  C              ::               output= at the sea-ice base to the ocean  C     iceMask     :: sea-ice fractional mask [0-1]
65  C     Tsf      :: surface (ice or snow) temperature  C     hIce        :: ice height [m]
66  C     qicen    :: ice enthalpy (J m-3)  C     hSnow       :: snow height [m]
67  C     Tice     :: internal ice temperatures  C     tFrz        :: sea-water freezing temperature [oC] (function of S)
68  C     sHeating :: surf heating left to melt snow or ice (= Atmos-conduction)  C     flxExSW     :: surf. heat flux (+=down) except SW, function of surf. temp Ts:
69  C     flxCnB   :: heat flux conducted through the ice to bottom surface  C                    0: Flx(Ts=0) ; 1: Flx(Ts=Ts^n) ; 2: d.Flx/dTs(Ts=Ts^n)
70  C     dTsf     :: surf. temp adjusment: Ts^n+1 - Ts^n  C---  Modified (input&output):
71  C     flxAtm   :: net flux of energy from the atmosphere [W/m2] (+=down)  C     flxSW       :: net Short-Wave flux (+=down) [W/m2]: input= at surface
72  C                without snow precip. (energy=0 for liquid water at 0.oC)  C                 ::                 output= below sea-ice, into the ocean
73  C     evpAtm   :: evaporation to the atmosphere (kg/m2/s) (>0 if evaporate)  C     tSrf        :: surface (ice or snow) temperature
74  C   i,j,bi,bj  :: indices of current grid point  C     qIc1        :: ice enthalpy (J/kg), 1rst level
75  C     myThid   :: Thread no. that called this routine.  C     qIc2        :: ice enthalpy (J/kg), 2nd level
76        LOGICAL useBlkFlx  C---  Output
77        _RL flxExcSw(0:2)  C     tIc1        :: temperature of ice layer 1 [oC]
78        _RL Tf  C     tIc2        :: temperature of ice layer 2 [oC]
79        _RL hi  C     dTsrf       :: surf. temp adjusment: Ts^n+1 - Ts^n
80        _RL hs  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        _RL flxSW  C     flxAtm      :: net flux of energy from the atmosphere [W/m2] (+=down)
83        _RL Tsf  C                    without snow precip. (energy=0 for liquid water at 0.oC)
84        _RL qicen(nlyr)  C     evpAtm      :: evaporation to the atmosphere [kg/m2/s] (>0 if evaporate)
85    C---  Input:
86        _RL Tice (nlyr)  C     myTime      :: current Time of simulation [s]
87        _RL sHeating  C     myIter      :: current Iteration number in simulation
88        _RL flxCnB  C     myThid      :: my Thread Id number
89        _RL dTsf        INTEGER bi,bj
90        _RL flxAtm        INTEGER iMin, iMax
91        _RL evpAtm        INTEGER jMin, jMax
92        INTEGER i,j, bi,bj        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        INTEGER myThid
114  CEOP  CEOP
115    
116  #ifdef ALLOW_THSICE  #ifdef ALLOW_THSICE
117    C     !LOCAL VARIABLES:
 C ADAPTED FROM:  
 C LANL CICE.v2.0.2  
 C-----------------------------------------------------------------------  
 C.. thermodynamics (vertical physics) based on M. Winton 3-layer model  
 C.. See Bitz, C. M. and W. H. Lipscomb, 1999:  "An energy-conserving  
 C..       thermodynamic sea ice model for climate study."  J. Geophys.  
 C..       Res., 104, 15669 - 15677.  
 C..     Winton, M., 1999:  "A reformulated three-layer sea ice model."    
 C..       Submitted to J. Atmos. Ocean. Technol.    
 C.. authors Elizabeth C. Hunke and William Lipscomb  
 C..         Fluid Dynamics Group, Los Alamos National Laboratory  
 C-----------------------------------------------------------------------  
 Cc****subroutine thermo_winton(n,fice,fsnow,dqice,dTsfc)  
 C.. Compute temperature change using Winton model with 2 ice layers, of  
 C.. which only the top layer has a variable heat capacity.  
   
118  C     == Local Variables ==  C     == Local Variables ==
119        INTEGER  k  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    C-    Define grid-point location where to print debugging values
174    #include "THSICE_DEBUG.h"
175    
176     1010 FORMAT(A,I3,3F11.6)
177     1020 FORMAT(A,1P4E14.6)
178    
179    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|
180    
181    #ifdef ALLOW_AUTODIFF_TAMC
182    c     act1 = bi - myBxLo(myThid)
183    c     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
184    c     act2 = bj - myByLo(myThid)
185    c     max2 = myByHi(myThid) - myByLo(myThid) + 1
186    c     act3 = myThid - 1
187    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
211    
212        _RL  frsnow        ! fractional snow cover        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        _RL  fswpen        ! SW penetrating beneath surface (W m-2)          ELSE
317        _RL  fswdn         ! SW absorbed at surface (W m-2)            iceFlag(i,j) = .FALSE.
318        _RL  fswint        ! SW absorbed in ice (W m-2)          ENDIF
319        _RL  fswocn        ! SW passed through ice to ocean (W m-2)         ENDDO
320          ENDDO
321        _RL  flxExceptSw   ! net surface heat flux, except short-wave (W/m2)        IF ( icount .gt. 0 ) THEN
322  C     evap           :: evaporation over snow/ice [kg/m2/s] (>0 if evaporate)         WRITE (standardMessageUnit,'(A,I5,A)')
323  C     dEvdT          :: derivative of evap. with respect to Tsf [kg/m2/s/K]       &      'THSICE_SOLVE4TEMP: there are ',icount,
324        _RL  evap, dEvdT       &      ' case(s) where hIce<hIceMin;'
325        _RL  flx0          ! net surf heat flux, from Atmos. to sea-ice (W m-2)         WRITE (standardMessageUnit,'(A,I3,A1,I3,A)')
326        _RL  fct           ! heat conducted to top surface       &      'THSICE_SOLVE4TEMP: the last one was at (',ii,',',jj,
327         &      ') with hIce = ', hIce(ii,jj)
328        _RL  df0dT         ! deriv of flx0 wrt Tsf (W m-2 deg-1)         STOP 'THSICE_SOLVE4TEMP: should not enter if hIce<hIceMin'
329          ENDIF
       _RL  k12, k32      ! thermal conductivity terms  
       _RL  a10, b10      ! coefficients in quadratic eqn for T1  
       _RL  a1, b1, c1    ! coefficients in quadratic eqn for T1  
 c     _RL  Tsf_start     ! old value of Tsf  
   
       _RL  dt            ! timestep  
   
       INTEGER iceornot  
       LOGICAL dBug  
   
  1010 FORMAT(A,I3,3F8.3)  
  1020 FORMAT(A,1P4E11.3)  
   
       dt  = thSIce_deltaT  
       dBug = .FALSE.  
 c     dBug = ( bi.EQ.3 .AND. i.EQ.15 .AND. j.EQ.11 )  
 c     dBug = ( bi.EQ.6 .AND. i.EQ.18 .AND. j.EQ.10 )  
       IF (dBug) WRITE(6,'(A,2I4,2I2)') 'ThSI_SOLVE4T: i,j=',i,j,bi,bj  
   
       if ( hi.lt.himin ) then  
 C If hi < himin, melt the ice.  
          STOP 'THSICE_SOLVE4TEMP: should not enter if hi<himin'  
       endif  
   
 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  
   
 C fractional snow cover  
       frsnow = 0. _d 0  
       if (hs .gt. 0. _d 0) frsnow = 1. _d 0  
   
 C Compute SW flux absorbed at surface and penetrating to layer 1.  
       fswpen  = flxSW * (1. _d 0 - frsnow) * i0  
       fswocn = fswpen * exp(-ksolar*hi)  
       fswint = fswpen - fswocn  
   
       fswdn = flxSW - 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 .or. Tice(2).gt.0. _d 0) then  
           write (6,*) 'BBerr Tice(1) > 0 = ',Tice(1)  
           write (6,*) 'BBerr Tice(2) > 0 = ',Tice(2)  
       endif  
       IF (dBug) WRITE(6,1010) 'ThSI_SOLVE4T: k, Ts, Tice=',0,Tsf,Tice  
   
 C Compute coefficients used in quadratic formula.  
   
       a10 = rhoi*cpice *hi/(2. _d 0*dt) +  
      &      k32 * (4. _d 0*dt*k32 + rhoi*cpice *hi)  
      &       / (6. _d 0*dt*k32 + rhoi*cpice *hi)  
       b10 = -hi*  
      &      (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.  
330    
331  C ----- begin iteration  -----  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|
       do 100 k = 1,nitMaxTsf  
332    
333  C Save temperatures at start of iteration.  #ifdef ALLOW_AUTODIFF_TAMC
334  c        Tsf_start = Tsf  CADJ STORE devdt    = comlev1_bibj,key=iicekey,byte=isbyte
335    CADJ STORE tsf      = comlev1_bibj,key=iicekey,byte=isbyte
336    #endif
337    
338          IF ( useBlkFlx ) THEN  C--   Get surface fluxes over melting surface
339  C Compute top surface flux.  #ifdef ALLOW_AUTODIFF_TAMC
340           if (hs.gt.3. _d -1) then        IF ( useBlkFlx ) THEN
341                iceornot=2  #else
342           else        IF ( useBlkFlx .AND. iterate4Tsf  ) THEN
343                iceornot=1  #endif
344           endif          DO j = jMin, jMax
345           CALL THSICE_GET_BULKF(           DO i = iMin, iMax
346       I                        iceornot, Tsf,             Tsf(i,j) = 0.
347       O                        flxExceptSw, df0dT, evap, dEvdT,           ENDDO
348       I                        i,j,bi,bj,myThid )          ENDDO
349           flx0 = fswdn + flxExceptSw          IF ( useEXF ) THEN
350          ELSE             CALL THSICE_GET_EXF(   bi, bj,
351           flx0 = fswdn + flxExcSw(1)       I                            iMin, iMax, jMin, jMax,
352           df0dT = flxExcSw(2)       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          ENDIF
365          IF ( dBug ) WRITE(6,1020) 'ThSI_SOLVE4T: flx0,df0dT,k12,D=',        ENDIF
366       &                                     flx0,df0dT,k12,k12-df0dT  
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
377            iterMax = nitMaxTsf
378          ELSE
379            iterMax = 1
380          ENDIF
381    
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 Compute new top layer and surface temperatures.  C ----- begin iteration  -----
390  C If Tsfc is computed to be > 0 C, fix Tsfc = 0 and recompute T1        DO k = 1,iterMax
391  C with different coefficients.  #ifndef ALLOW_AUTODIFF_TAMC
392           IF ( iterate4Tsf ) THEN
393           a1 = a10 - k12*df0dT / (k12-df0dT)         iterate4Tsf = .FALSE.
394           b1 = b10 - k12*(flx0-df0dT*Tsf) / (k12-df0dT)  #endif
395           Tice(1) = -(b1 + sqrt(b1*b1-4. _d 0*a1*c1))/(2. _d 0*a1)  
396           dTsf = (flx0 + k12*(Tice(1)-Tsf)) / (k12-df0dT)         IF ( useBlkFlx ) THEN
397           Tsf = Tsf + dTsf  C--   Compute top surface flux.
398           if (Tsf .gt. 0. _d 0) then           IF ( useEXF ) THEN
399              IF(dBug) WRITE(6,1010) 'ThSI_SOLVE4T: k,ts,t1,dTs=',             CALL THSICE_GET_EXF(   bi, bj,
400       &                                          k,Tsf,Tice(1),dTsf       I                            iMin, iMax, jMin, jMax,
401              a1 = a10 + k12       I                            iceFlag, hSnow, Tsf,
402              b1 = b10          ! note b1 = b10 - k12*Tf0       O                            flxTexSW, dFlxdT, evapT, dEvdT,
403              Tice(1) = (-b1 - sqrt(b1*b1-4. _d 0*a1*c1))/(2. _d 0*a1)       I                            myTime, myIter, myThid )
404              Tsf = 0. _d 0  C-    could add this "ifdef" to hide THSICE_GET_BULKF from TAF
405    #ifdef ALLOW_BULK_FORCE
406             ELSEIF ( useBulkForce ) THEN
407               CALL THSICE_GET_BULKF( bi, bj,
408         I                            iMin, iMax, jMin, jMax,
409         I                            iceFlag, hSnow, Tsf,
410         O                            flxTexSW, dFlxdT, evapT, dEvdT,
411         I                            myTime, myIter, myThid )
412    #endif /* ALLOW_BULK_FORCE */
413             ENDIF
414           ELSE
415             DO j = jMin, jMax
416              DO i = iMin, iMax
417               IF ( iceFlag(i,j) ) THEN
418                 flxTexSW(i,j) = flxExSW(i,j,1)
419                 dFlxdT(i,j) = flxExSW(i,j,2)
420               ENDIF
421              ENDDO
422             ENDDO
423           ENDIF
424    
425    #ifdef ALLOW_AUTODIFF_TAMC
426    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
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 )  
             flx0 = fswdn + flxExceptSw  
             dTsf = 0. _d 0  
468             ELSE             ELSE
469              flx0 = fswdn + 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              Tsf  = 0. _d 0            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 ) THEN  #ifdef ALLOW_DBUG_THSICE
485            if (abs(dTsf).lt.Terrmax) then            IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)
486              goto 150       &    'ThSI_SOLVE4T: k,ts,t1,dTs=', k,Tsf(i,j),tIc1(i,j),dTsrf(i,j)
487            elseif (k.eq.nitMaxTsf) then            IF ( useBlkFlx .AND. k.EQ.nitMaxTsf .AND. iceFlag(i,j) ) THEN
488              write (6,*) 'BB: thermw conv err ',i,j,bi,bj,dTsf             WRITE (6,'(A,4I4,I12,F15.9)')
489              write (6,*) 'BB: thermw conv err, iceheight ', hi       &       ' BB: not converge: i,j,it,hi=',i,j,bi,bj,myIter,hIce(i,j)
490              write (6,*) 'BB: thermw conv err: Tsf, flx0', Tsf,flx0             WRITE (6,*) 'BB: not converge: Tsf, dTsf=',
491              if (Tsf.lt.-70. _d 0) stop  !QQQQ fails       &          Tsf(i,j), dTsrf(i,j)
492            endif             WRITE (6,*) 'BB: not converge: flxNet,dFlxT=',
493           ELSE       &          flxNet, dFlxdT(i,j)
494              goto 150             IF (Tsf(i,j).LT.-70. _d 0) STOP
495           ENDIF            ENDIF
496    #endif
497    
498  100   continue  ! surface temperature iteration           ENDIF
499  150   continue          ENDDO
500           ENDDO
501    #ifndef ALLOW_AUTODIFF_TAMC
502           ENDIF
503    #endif
504          ENDDO
505  C ------ end iteration ------------  C ------ end iteration ------------
506    
507  c        if (Tsf.lt.-70. _d 0) then  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|
 cQQ        print*,'QQQQQ Tsf =',Tsf  
 cQQ        stop  !QQQQ fails  
 c        endif  
   
 C Compute new bottom layer temperature.  
   
       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  
508    
509          DO j = jMin, jMax
510  C Compute final flux values at surfaces.         DO i = iMin, iMax
511            IF ( iceMask(i,j).GT.0. _d 0 ) THEN
512        fct    = k12*(Tsf-Tice(1))  
513        flxCnB = 4. _d 0*kice *(Tice(2)-Tf)/hi  C--   Compute new bottom layer temperature.
514        flx0   = flx0 + df0dT*dTsf            k32       = 2. _d 0*kIce  / hIce(i,j)
515        IF ( useBlkFlx ) THEN            tIc2(i,j) = ( 2. _d 0*dt*k32*(tIc1(i,j)+2. _d 0*tFrz(i,j))
516          evpAtm = evap       &                 + rhoi*cpIce*hIce(i,j)*tIc2(i,j))
517  C--   needs to update also Evap (Tsf changes) since Latent heat has been updated       &               /(6. _d 0*dt*k32 + rhoi*cpIce*hIce(i,j))
518          evpAtm = evap + dEvdT*dTsf  #ifdef ALLOW_DBUG_THSICE
519  C-    energy flux to Atmos: use net short-wave flux at surf. and            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)
535                evpAtm(i,j) = 0.
536              ENDIF
537    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            flxAtm(i,j) = flxAtm(i,j) + flxTexSW(i,j)
540       &                 + evpAtm*Lfresh       &                + dFlxdT(i,j)*dTsrf(i,j) + evpAtm(i,j)*Lfresh
541        ELSE  C-    excess of energy @ surface (used for surface melting):
542          flxAtm = 0.            sHeat(i,j) = flxNet - fct
543          evpAtm = 0.  
544        ENDIF  #ifdef ALLOW_DBUG_THSICE
545        sHeating = flx0 - fct            IF ( dBug(i,j,bi,bj) ) WRITE(6,1020)
546         &     'ThSI_SOLVE4T: flxNet,fct,Dif,flxCnB=',
547  C-    SW flux at sea-ice base left to the ocean       &                    flxNet,fct,flxNet-fct,flxCnB(i,j)
548        flxSW = fswocn  #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        IF (dBug) WRITE(6,1020) 'ThSI_SOLVE4T: flx0,fct,Dif,flxCnB=',          ELSE
578       &    flx0,fct,flx0-fct,flxCnB  C--     ice-free grid point:
579    c         tIc1  (i,j) = 0. _d 0
580  C Compute new enthalpy for each layer.  c         tIc2  (i,j) = 0. _d 0
581              dTsrf (i,j) = 0. _d 0
582        qicen(1) = -cpwater*Tmlt1 + cpice *(Tmlt1-Tice(1)) +  c         sHeat (i,j) = 0. _d 0
583       &              Lfresh*(1. _d 0-Tmlt1/Tice(1))  c         flxCnB(i,j) = 0. _d 0
584        qicen(2) = -cpice *Tice(2) + Lfresh  c         flxAtm(i,j) = 0. _d 0
585    c         evpAtm(i,j) = 0. _d 0
 C Make sure internal ice temperatures do not exceed Tmlt.  
 C (This should not happen for reasonable values of i0.)  
   
       if (Tice(1) .ge. Tmlt1) then  
         write (6,'(A,2I4,2I3,1P2E14.6)')  
      &   'BBerr - Bug: IceT(1) > Tmlt',i,j,bi,bj,Tice(1),Tmlt1  
       endif  
       if (Tice(2) .ge. 0. _d 0) then  
        write (6,'(A,2I4,2I3,1P2E14.6)')  
      &   'BBerr - Bug: IceT(2) > 0',i,j,bi,bj,Tice(2)  
       endif  
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.5  
changed lines
  Added in v.1.24
  ViewVC Help
Powered by ViewVC 1.1.22