/[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.1 by jmc, Wed Apr 7 23:40:34 2004 UTC revision 1.14 by heimbach, Mon Apr 16 22:38:24 2007 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, siLo, siHi, sjLo, sjHi,
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."  J. Geophys.
32    C..       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"
47  #include "THSICE_SIZE.h"  #include "THSICE_SIZE.h"
48  #include "THSICE_PARAMS.h"  #include "THSICE_PARAMS.h"
49    #ifdef ALLOW_AUTODIFF_TAMC
50    # include "SIZE.h"
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     siLo,siHi   :: size of input/output array: 1rst dim. lower,higher bounds
58  C     flxExcSw :: surf. heat flux (+=down) except SW, function of surf. temp Ts:  C     sjLo,sjHi   :: size of input/output array: 2nd  dim. lower,higher bounds
59  C                0: Flx(Ts=0) ; 1: Flx(Ts=Ts^n) ; 2: d.Flx/dTs(Ts=Ts^n)  C     bi,bj       :: tile indices
60  C     Tf       :: freezing temperature (oC) of local sea-water  C     iMin,iMax   :: computation domain: 1rst index range
61  C     hi       :: ice height  C     jMin,jMax   :: computation domain: 2nd  index range
62  C     hs       :: snow height  C     dBugFlag    :: allow to print debugging stuff (e.g. on 1 grid point).
63  C     flxSW    :: net Short-Wave flux (+=down) [W/m2]: input= at surface  C     useBulkForce:: use surf. fluxes from bulk-forcing external S/R
64  C              ::               output= at the sea-ice base to the ocean  C     useEXF      :: use surf. fluxes from exf          external S/R
65  C     Tsf      :: surface (ice or snow) temperature  C---  Input:
66  C     qicen    :: ice enthalpy (J m-3)  C         iceMask :: sea-ice fractional mask [0-1]
67  C     Tice     :: internal ice temperatures  C  hIce    (hi)   :: ice height [m]
68  C     sHeating :: surf heating left to melt snow or ice (= Atmos-conduction)  C  hSnow   (hs)   :: snow height [m]
69  C     flxCnB   :: heat flux conducted through the ice to bottom surface  C  tFrz    (Tf)   :: sea-water freezing temperature [oC] (function of S)
70  C     dTsf     :: surf. temp adjusment: Ts^n+1 - Ts^n  C  flxExSW  (=)   :: surf. heat flux (+=down) except SW, function of surf. temp Ts:
71  C     flxAtm   :: net flux of energy from the atmosphere [W/m2] (+=down)  C                    0: Flx(Ts=0) ; 1: Flx(Ts=Ts^n) ; 2: d.Flx/dTs(Ts=Ts^n)
72  C                without snow precip. (energy=0 for liquid water at 0.oC)  C---  Modified (input&output):
73  C     evpAtm   :: evaporation to the atmosphere (kg/m2/s) (>0 if evaporate)  C  flxSW  (netSW) :: net Short-Wave flux (+=down) [W/m2]: input= at surface
74  C   i,j,bi,bj  :: indices of current grid point  C           (=)   ::                 output= below sea-ice, into the ocean
75  C     myThid   :: Thread no. that called this routine.  C  tSrf    (Tsf)  :: surface (ice or snow) temperature
76        LOGICAL useBlkFlx  C  qIc1   (qicen) :: ice enthalpy (J/kg), 1rst level
77        _RL flxExcSw(0:2)  C  qIc2   (qicen) :: ice enthalpy (J/kg), 2nd level
78    C---  Output
79    C  tIc1    (Tice) :: temperature of ice layer 1 [oC]
80    C  tIc2    (Tice) :: temperature of ice layer 2 [oC]
81    C  dTsrf   (dTsf) :: surf. temp adjusment: Ts^n+1 - Ts^n
82    C  sHeat(sHeating):: surf heating flux left to melt snow or ice (= Atmos-conduction)
83    C  flxCnB   (=)   :: heat flux conducted through the ice to bottom surface
84    C  flxAtm   (=)   :: net flux of energy from the atmosphere [W/m2] (+=down)
85    C                    without snow precip. (energy=0 for liquid water at 0.oC)
86    C  evpAtm   (=)   :: evaporation to the atmosphere [kg/m2/s] (>0 if evaporate)
87    C---  Input:
88    C     myTime      :: current Time of simulation [s]
89    C     myIter      :: current Iteration number in simulation
90    C     myThid      :: my Thread Id number
91          INTEGER siLo, siHi, sjLo, sjHi
92          INTEGER bi,bj
93          INTEGER iMin, iMax
94          INTEGER jMin, jMax
95          LOGICAL dBugFlag
96          LOGICAL useBulkForce
97          LOGICAL useEXF
98          _RL iceMask(siLo:siHi,sjLo:sjHi)
99          _RL hIce   (siLo:siHi,sjLo:sjHi)
100          _RL hSnow  (siLo:siHi,sjLo:sjHi)
101          _RL tFrz   (siLo:siHi,sjLo:sjHi)
102          _RL flxExSW(iMin:iMax,jMin:jMax,0:2)
103          _RL flxSW  (siLo:siHi,sjLo:sjHi)
104          _RL tSrf   (siLo:siHi,sjLo:sjHi)
105          _RL qIc1   (siLo:siHi,sjLo:sjHi)
106          _RL qIc2   (siLo:siHi,sjLo:sjHi)
107          _RL tIc1   (siLo:siHi,sjLo:sjHi)
108          _RL tIc2   (siLo:siHi,sjLo:sjHi)
109    c     _RL dTsrf  (siLo:siHi,sjLo:sjHi)
110          _RL dTsrf  (iMin:iMax,jMin:jMax)
111          _RL sHeat  (siLo:siHi,sjLo:sjHi)
112          _RL flxCnB (siLo:siHi,sjLo:sjHi)
113          _RL flxAtm (siLo:siHi,sjLo:sjHi)
114          _RL evpAtm (siLo:siHi,sjLo:sjHi)
115          _RL  myTime
116          INTEGER myIter
117          INTEGER myThid
118    CEOP
119    
120    #ifdef ALLOW_THSICE
121    C     !LOCAL VARIABLES:
122    C---  local copy of input/output argument list variables (see description above)
123    c     _RL flxExcSw(0:2)
124        _RL Tf        _RL Tf
125        _RL hi        _RL hi
126        _RL hs        _RL hs
127          _RL netSW
       _RL flxSW  
128        _RL Tsf        _RL Tsf
129        _RL qicen(nlyr)        _RL qicen(nlyr)
   
130        _RL Tice (nlyr)        _RL Tice (nlyr)
131        _RL sHeating  c     _RL sHeating
132        _RL flxCnB  c     _RL flxCnB
133        _RL dTsf        _RL dTsf
134        _RL flxAtm  c     _RL flxAtm
135        _RL evpAtm  c     _RL evpAtm
       INTEGER i,j, bi,bj  
       INTEGER myThid  
 CEOP  
   
 #ifdef ALLOW_THSICE  
   
 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.  
136    
137  C     == Local Variables ==  C     == Local Variables ==
138        INTEGER  k        INTEGER i,j
139          INTEGER k, iterMax
140    
141        _RL  frsnow        ! fractional snow cover        _RL  frsnow        ! fractional snow cover
   
142        _RL  fswpen        ! SW penetrating beneath surface (W m-2)        _RL  fswpen        ! SW penetrating beneath surface (W m-2)
143        _RL  fswdn         ! SW absorbed at surface (W m-2)        _RL  fswdn         ! SW absorbed at surface (W m-2)
144        _RL  fswint        ! SW absorbed in ice (W m-2)        _RL  fswint        ! SW absorbed in ice (W m-2)
145        _RL  fswocn        ! SW passed through ice to ocean (W m-2)        _RL  fswocn        ! SW passed through ice to ocean (W m-2)
   
146        _RL  flxExceptSw   ! net surface heat flux, except short-wave (W/m2)        _RL  flxExceptSw   ! net surface heat flux, except short-wave (W/m2)
147  C     evap           :: evaporation over snow/ice [kg/m2/s] (>0 if evaporate)  C     evap           :: evaporation over snow/ice [kg/m2/s] (>0 if evaporate)
148  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]
149        _RL  evap, dEvdT        _RL  evap, dEvdT
150        _RL  flx0          ! net surf heat flux, from Atmos. to sea-ice (W m-2)        _RL  flx0          ! net surf heat flux, from Atmos. to sea-ice (W m-2)
151        _RL  fct           ! heat conducted to top surface        _RL  fct           ! heat conducted to top surface
   
152        _RL  df0dT         ! deriv of flx0 wrt Tsf (W m-2 deg-1)        _RL  df0dT         ! deriv of flx0 wrt Tsf (W m-2 deg-1)
   
153        _RL  k12, k32      ! thermal conductivity terms        _RL  k12, k32      ! thermal conductivity terms
154        _RL  a10, b10      ! coefficients in quadratic eqn for T1        _RL  a10, b10      ! coefficients in quadratic eqn for T1
155        _RL  a1, b1, c1    ! coefficients in quadratic eqn for T1        _RL  a1, b1, c1    ! coefficients in quadratic eqn for T1
156  c     _RL  Tsf_start     ! old value of Tsf  c     _RL  Tsf_start     ! old value of Tsf
   
157        _RL  dt            ! timestep        _RL  dt            ! timestep
   
158        INTEGER iceornot        INTEGER iceornot
159        LOGICAL dBug        LOGICAL useBlkFlx
160    
161    C-    define grid-point location where to print debugging values
162    #include "THSICE_DEBUG.h"
163    
164   1010 FORMAT(A,I3,3F8.3)   1010 FORMAT(A,I3,3F11.6)
165   1020 FORMAT(A,1P4E11.3)   1020 FORMAT(A,1P4E14.6)
166    
167        dt  = thSIce_deltaT  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
       dBug = .FALSE.  
 c     dBug = ( bi.EQ.3 .AND. i.EQ.13 .AND. j.EQ.13 )  
 c     dBug = ( bi.EQ.6 .AND. i.EQ.18 .AND. j.EQ.10 )  
       IF ( dBug ) WRITE(6,'(A,2I4,2I2)') 'ThSI_THERM: i,j=',i,j,bi,bj  
168    
169        if ( hi.lt.himin ) then  #ifdef ALLOW_AUTODIFF_TAMC
170              act1 = bi - myBxLo(myThid)
171              max1 = myBxHi(myThid) - myBxLo(myThid) + 1
172              act2 = bj - myByLo(myThid)
173              max2 = myByHi(myThid) - myByLo(myThid) + 1
174              act3 = myThid - 1
175              max3 = nTx*nTy
176              act4 = ikey_dynamics - 1
177    #endif /* ALLOW_AUTODIFF_TAMC */
178    
179          useBlkFlx = useEXF .OR. useBulkForce
180    
181          dt  = thSIce_dtTemp
182          DO j = jMin, jMax
183           DO i = iMin, iMax
184    #ifdef ALLOW_AUTODIFF_TAMC
185              ikey_1 = i
186         &         + sNx*(j-1)
187         &         + sNx*sNy*act1
188         &         + sNx*sNy*max1*act2
189         &         + sNx*sNy*max1*max2*act3
190         &         + sNx*sNy*max1*max2*max3*act4
191    #endif /* ALLOW_AUTODIFF_TAMC */
192    C--
193    #ifdef ALLOW_AUTODIFF_TAMC
194    CADJ STORE  df0dt      = comlev1_thsice_1, key=ikey_1
195    CADJ STORE  flxsw(i,j) = comlev1_thsice_1, key=ikey_1
196    CADJ STORE  qic1(i,j)  = comlev1_thsice_1, key=ikey_1
197    CADJ STORE  qic2(i,j)  = comlev1_thsice_1, key=ikey_1
198    CADJ STORE  tsrf(i,j)  = comlev1_thsice_1, key=ikey_1
199    #endif
200            IF ( iceMask(i,j).GT.0. _d 0) THEN
201              hi      = hIce(i,j)
202              hs      = hSnow(i,j)
203              Tf      = tFrz(i,j)
204              netSW   = flxSW(i,j)
205              Tsf     = tSrf(i,j)
206              qicen(1)= qIc1(i,j)
207              qicen(2)= qIc2(i,j)
208    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
209    #ifdef ALLOW_DBUG_THSICE
210              IF ( dBug(i,j,bi,bj) ) WRITE(6,'(A,2I4,2I2)')
211         &         'ThSI_SOLVE4T: i,j=',i,j,bi,bj
212    #endif
213          IF ( hi.LT.himin ) THEN
214  C If hi < himin, melt the ice.  C If hi < himin, melt the ice.
215           STOP 'THSICE_THERM: should not enter if hi<himin'           STOP 'THSICE_SOLVE4TEMP: should not enter if hi<himin'
216        endif        ENDIF
217    
218  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
219    
220  C fractional snow cover  C fractional snow cover
221        frsnow = 0. _d 0        frsnow = 0. _d 0
222        if (hs .gt. 0. _d 0) frsnow = 1. _d 0        IF (hs .GT. 0. _d 0) frsnow = 1. _d 0
223    
224  C Compute SW flux absorbed at surface and penetrating to layer 1.  C Compute SW flux absorbed at surface and penetrating to layer 1.
225        fswpen  = flxSW * (1. _d 0 - frsnow) * i0        fswpen  = netSW * (1. _d 0 - frsnow) * i0
226        fswocn = fswpen * exp(-ksolar*hi)        fswocn = fswpen * exp(-ksolar*hi)
227        fswint = fswpen - fswocn        fswint = fswpen - fswocn
228    
229        fswdn = flxSW - fswpen        fswdn = netSW - fswpen
230    
231  C Compute conductivity terms at layer interfaces.  C Compute conductivity terms at layer interfaces.
232    
# Line 151  C compute ice temperatures Line 237  C compute ice temperatures
237        a1 = cpice        a1 = cpice
238        b1 = qicen(1) + (cpwater-cpice )*Tmlt1 - Lfresh        b1 = qicen(1) + (cpwater-cpice )*Tmlt1 - Lfresh
239        c1 = Lfresh * Tmlt1        c1 = Lfresh * Tmlt1
240        Tice(1) = 0.5 _d 0 *(-b1 - sqrt(b1*b1-4. _d 0*a1*c1))/a1        Tice(1) = 0.5 _d 0 *(-b1 - SQRT(b1*b1-4. _d 0*a1*c1))/a1
241        Tice(2) = (Lfresh-qicen(2)) / cpice        Tice(2) = (Lfresh-qicen(2)) / cpice
242    
243        if (Tice(1).gt.0. _d 0 .or. Tice(2).gt.0. _d 0) then        IF (Tice(1).GT.0. _d 0 ) THEN
244            write (6,*) 'BBerr Tice(1) > 0 = ',Tice(1)         WRITE (standardMessageUnit,'(A,I12,1PE14.6)')
245            write (6,*) 'BBerr Tice(2) > 0 = ',Tice(2)       &       ' BBerr: Tice(1) > 0 ; it=', myIter, qicen(1)
246        endif         WRITE (standardMessageUnit,'(A,4I5,2F11.4)')
247        IF ( dBug ) WRITE(6,1010) 'ThSI_THERM: k,Tice=',0,Tsf,Tice       &      ' BBerr: i,j,bi,bj,Tice = ',i,j,bi,bj,Tice
248  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|        ENDIF
249          IF ( Tice(2).GT.0. _d 0) THEN
250           WRITE (standardMessageUnit,'(A,I12,1PE14.6)')
251         &       ' BBerr: Tice(2) > 0 ; it=', myIter, qicen(2)
252           WRITE (standardMessageUnit,'(A,4I5,2F11.4)')
253         &      ' BBerr: i,j,bi,bj,Tice = ',i,j,bi,bj,Tice
254          ENDIF
255    #ifdef ALLOW_DBUG_THSICE
256          IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)
257         &  'ThSI_SOLVE4T: k, Ts, Tice=',0,Tsf,Tice
258    #endif
259    
260  C Compute coefficients used in quadratic formula.  C Compute coefficients used in quadratic formula.
261    
# Line 173  C Compute coefficients used in quadratic Line 269  C Compute coefficients used in quadratic
269       &       / (6. _d 0*dt*k32 + rhoi*cpice *hi) - fswint       &       / (6. _d 0*dt*k32 + rhoi*cpice *hi) - fswint
270        c1 = rhoi*Lfresh*hi*Tmlt1 / (2. _d 0*dt)        c1 = rhoi*Lfresh*hi*Tmlt1 / (2. _d 0*dt)
271    
272    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
273  C Compute new surface and internal temperatures; iterate until  C Compute new surface and internal temperatures; iterate until
274  C Tsfc converges.  C Tsfc converges.
275    
276          IF ( useBlkFlx ) THEN
277            iterMax = nitMaxTsf
278          ELSE
279            iterMax = 1
280          ENDIF
281          dTsf = Terrmax
282    
283  C ----- begin iteration  -----  C ----- begin iteration  -----
284        do 100 k = 1,nitMaxTsf        DO k = 1,iterMax
285    
286    #ifdef ALLOW_AUTODIFF_TAMC
287           ikey_3 = (ikey_1-1)*MaxTsf + k
288    #endif
289    
290    #ifdef ALLOW_AUTODIFF_TAMC
291    CADJ STORE dtsf = comlev1_thsice_3, key=ikey_3
292    #endif
293           IF ( ABS(dTsf).GE.Terrmax ) THEN
294    
295  C Save temperatures at start of iteration.  C Save temperatures at start of iteration.
296  c        Tsf_start = Tsf  c        Tsf_start = Tsf
297    
298          IF ( useBlkFlx ) THEN          IF ( useBlkFlx ) THEN
299  C Compute top surface flux.  C Compute top surface flux.
300           if (hs.gt.3. _d -1) then           IF (hs.GT.3. _d -1) THEN
301                iceornot=2                iceornot=2
302           else           ELSE
303                iceornot=1                iceornot=1
304           endif           ENDIF
305           CALL THSICE_GET_BULKF(           IF ( useBulkForce ) THEN
306       I                        iceornot, Tsf,            CALL THSICE_GET_BULKF(
307       O                        flxExceptSw, df0dT, evap, dEvdT,       I                         iceornot, Tsf,
308       I                        i,j,bi,bj,myThid )       O                         flxExceptSw, df0dT, evap, dEvdT,
309           flx0 = fswdn + flxExceptSw       I                         i,j,bi,bj,myThid )
310             ELSEIF ( useEXF ) THEN
311              CALL THSICE_GET_EXF  (
312         I                         iceornot, Tsf,
313         O                         flxExceptSw, df0dT, evap, dEvdT,
314         I                         i,j,bi,bj,myThid )
315             ENDIF
316          ELSE          ELSE
317           flx0 = fswdn + flxExcSw(1)           flxExceptSw = flxExSW(i,j,1)
318           df0dT = flxExcSw(2)           df0dT = flxExSW(i,j,2)
319          ENDIF          ENDIF
320            flx0 = fswdn + flxExceptSw
321    #ifdef ALLOW_DBUG_THSICE
322          IF ( dBug(i,j,bi,bj) ) WRITE(6,1020)
323         &  'ThSI_SOLVE4T: flx0,df0dT,k12,D=', flx0,df0dT,k12,k12-df0dT
324    #endif
325    
326  C Compute new top layer and surface temperatures.  C Compute new top layer and surface temperatures.
327  C If Tsfc is computed to be > 0 C, fix Tsfc = 0 and recompute T1  C If Tsfc is computed to be > 0 C, fix Tsfc = 0 and recompute T1
328  C with different coefficients.  C with different coefficients.
329    
330    #ifdef ALLOW_AUTODIFF_TAMC
331    CADJ STORE tsf  = comlev1_thsice_3, key=ikey_3
332    #endif
333           a1 = a10 - k12*df0dT / (k12-df0dT)           a1 = a10 - k12*df0dT / (k12-df0dT)
334           b1 = b10 - k12*(flx0-df0dT*Tsf) / (k12-df0dT)           b1 = b10 - k12*(flx0-df0dT*Tsf) / (k12-df0dT)
335           Tice(1) = -(b1 + sqrt(b1*b1-4. _d 0*a1*c1))/(2. _d 0*a1)           Tice(1) = -(b1 + SQRT(b1*b1-4. _d 0*a1*c1))/(2. _d 0*a1)
336           dTsf = (flx0 + k12*(Tice(1)-Tsf)) / (k12-df0dT)           dTsf = (flx0 + k12*(Tice(1)-Tsf)) / (k12-df0dT)
337           Tsf = Tsf + dTsf           Tsf = Tsf + dTsf
338           if (Tsf .gt. 0. _d 0) then           IF (Tsf .GT. 0. _d 0) THEN
339              IF(dBug) WRITE(6,1010) 'ThSI_THERM: k,tice=',  #ifdef ALLOW_DBUG_THSICE
340       &                                          k,Tsf,Tice(1),dTsf              IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)
341         &        'ThSI_SOLVE4T: k,ts,t1,dTs=', k,Tsf,Tice(1),dTsf
342    #endif
343              a1 = a10 + k12              a1 = a10 + k12
344              b1 = b10          ! note b1 = b10 - k12*Tf0              b1 = b10          ! note b1 = b10 - k12*Tf0
345              Tice(1) = (-b1 - sqrt(b1*b1-4. _d 0*a1*c1))/(2. _d 0*a1)              Tice(1) = (-b1 - SQRT(b1*b1-4. _d 0*a1*c1))/(2. _d 0*a1)
346              Tsf = 0. _d 0              Tsf = 0. _d 0
347             IF ( useBlkFlx ) THEN             IF ( useBlkFlx ) THEN
348              if (hs.gt.3. _d -1) then              IF (hs.GT.3. _d -1) THEN
349                   iceornot=2                   iceornot=2
350              else              ELSE
351                   iceornot=1                   iceornot=1
352              endif              ENDIF
353              CALL THSICE_GET_BULKF(              IF ( useBulkForce ) THEN
354       I                        iceornot, Tsf,               CALL THSICE_GET_BULKF(
355       O                        flxExceptSw, df0dT, evap, dEvdT,       I                            iceornot, Tsf,
356       I                        i,j,bi,bj,myThid )       O                            flxExceptSw, df0dT, evap, dEvdT,
357              flx0 = fswdn + flxExceptSw       I                            i,j,bi,bj,myThid )
358                ELSEIF ( useEXF ) THEN
359                 CALL THSICE_GET_EXF  (
360         I                            iceornot, Tsf,
361         O                            flxExceptSw, df0dT, evap, dEvdT,
362         I                            i,j,bi,bj,myThid )
363                ENDIF
364              dTsf = 0. _d 0              dTsf = 0. _d 0
365             ELSE             ELSE
366              flx0 = fswdn + flxExcSw(0)              flxExceptSw = flxExSW(i,j,0)
367              dTsf = 1000.              dTsf = 1000.
368              df0dT = 0.              df0dT = 0.
369             ENDIF             ENDIF
370              Tsf  = 0. _d 0             flx0 = fswdn + flxExceptSw
371           endif           ENDIF
372    
373  C Check for convergence.  If no convergence, then repeat.  C Check for convergence.  If no convergence, then repeat.
374  C  C
375  C Convergence test: Make sure Tsfc has converged, within prescribed error.    C Convergence test: Make sure Tsfc has converged, within prescribed error.
376  C (Energy conservation is guaranteed within machine roundoff, even  C (Energy conservation is guaranteed within machine roundoff, even
377  C if Tsfc has not converged.)  C if Tsfc has not converged.)
378  C If no convergence, then repeat.  C If no convergence, then repeat.
379    
380           IF ( dBug ) WRITE(6,1010) 'ThSI_THERM: k,tice=',  #ifdef ALLOW_DBUG_THSICE
381       &                                            k,Tsf,Tice(1),dTsf           IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)
382           IF ( useBlkFlx ) THEN       &     'ThSI_SOLVE4T: k,ts,t1,dTs=', k,Tsf,Tice(1),dTsf
383            if (abs(dTsf).lt.Terrmax) then  #endif
384              goto 150           IF ( useBlkFlx .AND. k.EQ.nitMaxTsf
385            elseif (k.eq.nitMaxTsf) then       &                  .AND. ABS(dTsf).GE.Terrmax ) THEN
386              write (6,*) 'BB: thermw conv err ',i,j,bi,bj,dTsf              WRITE (6,'(A,4I4,I12,F15.9)')
387              write (6,*) 'BB: thermw conv err, iceheight ', hi       &                 ' BB: not converge: i,j,it,hi=',i,j,bi,bj,
388              write (6,*) 'BB: thermw conv err: Tsf, flx0', Tsf,flx0       &                                                    myIter,hi
389              if (Tsf.lt.-70. _d 0) stop  !QQQQ fails              WRITE (6,*) 'BB: not converge: Tsf, dTsf=', Tsf,dTsf
390            endif              WRITE (6,*) 'BB: not converge: flx0,dfdT=',flx0,df0dT
391           ELSE              IF (Tsf.LT.-70. _d 0) STOP
             goto 150  
392           ENDIF           ENDIF
393    
394  100   continue  ! surface temperature iteration         ENDIF
395  150   continue        ENDDO
396  C ------ end iteration ------------  C ------ end iteration ------------
397    
 c        if (Tsf.lt.-70. _d 0) then  
 cQQ        print*,'QQQQQ Tsf =',Tsf  
 cQQ        stop  !QQQQ fails  
 c        endif  
   
398  C Compute new bottom layer temperature.  C Compute new bottom layer temperature.
399    
400        Tice(2) = (2. _d 0*dt*k32*(Tice(1)+2. _d 0*Tf)        Tice(2) = (2. _d 0*dt*k32*(Tice(1)+2. _d 0*Tf)
401       &        + rhoi*cpice *hi*Tice(2))       &        + rhoi*cpice *hi*Tice(2))
402       &         /(6. _d 0*dt*k32 + rhoi*cpice *hi)       &         /(6. _d 0*dt*k32 + rhoi*cpice *hi)
403        IF ( dBug ) WRITE(6,1010) 'ThSI_THERM: k,Tice=',k,Tsf,Tice  #ifdef ALLOW_DBUG_THSICE
404          IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)
405         &  'ThSI_SOLVE4T: k, Ts, Tice=',k,Tsf,Tice
406    #endif
407    
408  C Compute final flux values at surfaces.  C Compute final flux values at surfaces.
409    
410        fct    = k12*(Tsf-Tice(1))        fct    = k12*(Tsf-Tice(1))
411        flxCnB = 4. _d 0*kice *(Tice(2)-Tf)/hi        flxCnB(i,j) = 4. _d 0*kice *(Tice(2)-Tf)/hi
412        flx0   = flx0 + df0dT*dTsf        flx0   = flx0 + df0dT*dTsf
413        IF ( useBlkFlx ) THEN        IF ( useBlkFlx ) THEN
         evpAtm = evap  
414  C--   needs to update also Evap (Tsf changes) since Latent heat has been updated  C--   needs to update also Evap (Tsf changes) since Latent heat has been updated
415  C     Note: this fix affects the results and is postponed for now.          evpAtm(i,j) = evap + dEvdT*dTsf
 c       evpAtm = evap + dEvdT*dTsf  
 C-    energy flux to Atmos: use net short-wave flux at surf. and  
 C     use latent heat = Lvap (energy=0 for liq. water at 0.oC)  
         flxAtm = flxSW + flxExceptSw + df0dT*dTsf  
      &                 + evpAtm*Lfresh  
416        ELSE        ELSE
417          flxAtm = 0.  C- WARNING: Evap & +Evap*Lfresh are missing ! (but only affects Diagnostics)
418          evpAtm = 0.          evpAtm(i,j) = 0.
419        ENDIF        ENDIF
420        sHeating = flx0 - fct  C-    energy flux to Atmos: use net short-wave flux at surf. and
421    C     use latent heat = Lvap (energy=0 for liq. water at 0.oC)
422          flxAtm(i,j) = netSW + flxExceptSw
423         &                    + df0dT*dTsf + evpAtm(i,j)*Lfresh
424    C-    excess of energy @ surface (used for surface melting):
425          sHeat(i,j) = flx0 - fct
426    
427  C-    SW flux at sea-ice base left to the ocean  C-    SW flux at sea-ice base left to the ocean
428        flxSW = fswocn        flxSW(i,j) = fswocn
429    
430        IF (dBug) WRITE(6,1020) 'ThSI_THERM: flx0,fct,Dif,flxCnB=',  #ifdef ALLOW_DBUG_THSICE
431       &    flx0,fct,flx0-fct,flxCnB        IF ( dBug(i,j,bi,bj) ) WRITE(6,1020)
432         &  'ThSI_SOLVE4T: flx0,fct,Dif,flxCnB=',
433         &                 flx0,fct,flx0-fct,flxCnB(i,j)
434    #endif
435    
436  C Compute new enthalpy for each layer.  C Compute new enthalpy for each layer.
437    
438        qicen(1) = -cpwater*Tmlt1 + cpice *(Tmlt1-Tice(1)) +        qicen(1) = -cpwater*Tmlt1 + cpice *(Tmlt1-Tice(1))
439       &              Lfresh*(1. _d 0-Tmlt1/Tice(1))       &            + Lfresh*(1. _d 0-Tmlt1/Tice(1))
440        qicen(2) = -cpice *Tice(2) + Lfresh        qicen(2) = -cpice *Tice(2) + Lfresh
441    
442  C Make sure internal ice temperatures do not exceed Tmlt.  C Make sure internal ice temperatures do not exceed Tmlt.
443  C (This should not happen for reasonable values of i0.)  C (This should not happen for reasonable values of i0.)
444    
445        if (Tice(1) .ge. Tmlt1) then        IF (Tice(1) .GE. Tmlt1) THEN
446          write (6,'(A,2I4,2I3,1P2E14.6)')          WRITE (6,'(A,2I4,2I3,1P2E14.6)')
447       &   'BBerr - Bug: IceT(1) > Tmlt',i,j,bi,bj,Tice(1),Tmlt1       &   ' BBerr - Bug: IceT(1) > Tmlt',i,j,bi,bj,Tice(1),Tmlt1
448        endif        ENDIF
449        if (Tice(2) .ge. 0. _d 0) then        IF (Tice(2) .GE. 0. _d 0) THEN
450         write (6,'(A,2I4,2I3,1P2E14.6)')         WRITE (6,'(A,2I4,2I3,1P2E14.6)')
451       &   'BBerr - Bug: IceT(2) > 0',i,j,bi,bj,Tice(2)       &   ' BBerr - Bug: IceT(2) > 0',i,j,bi,bj,Tice(2)
452        endif        ENDIF
453    
454    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
455    C--    Update Sea-Ice state :
456              tSrf(i,j)  = Tsf
457              tIc1(i,j)  = Tice(1)
458              tic2(i,j)  = Tice(2)
459              qIc1(i,j)  = qicen(1)
460              qIc2(i,j)  = qicen(2)
461    c         dTsrf(i,j) = dTsf
462              IF ( .NOT.useBlkFlx ) dTsrf(i,j) = dTsf
463    c         sHeat(i,j) = sHeating
464    c         flxCnB(i,j)= flxCnB
465    c         flxAtm(i,j)= flxAtm
466    c         evpAtm(i,j)= evpAtm
467    #ifdef ALLOW_DBUG_THSICE
468              IF ( dBug(i,j,bi,bj) ) THEN
469               WRITE(6,1020) 'ThSI_SOLV_4T: Tsf, Tice(1,2), dTsurf=',
470         &                              Tsf, Tice, dTsf
471               WRITE(6,1020) 'ThSI_SOLV_4T: sHeat, flxCndBt, Qice =',
472         &           sHeat(i,j), flxCnB(i,j), qicen
473               WRITE(6,1020) 'ThSI_SOLV_4T: flxA, evpA, fxSW_bf,af=',
474         &           flxAtm(i,j), evpAtm(i,j), netSW, flxSW(i,j)
475              ENDIF
476    #endif
477            ELSE
478              IF ( .NOT.useBlkFlx ) dTsrf(i,j) = 0. _d 0
479            ENDIF
480           ENDDO
481          ENDDO
482  #endif  /* ALLOW_THSICE */  #endif  /* ALLOW_THSICE */
483    
484  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.14
  ViewVC Help
Powered by ViewVC 1.1.22