/[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.16 by jmc, Sun Apr 29 23:48:44 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    
# Line 27  C     == Global variables === Line 46  C     == Global variables ===
46  #include "EEPARAMS.h"  #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   1010 FORMAT(A,I3,3F8.3)  C-    define grid-point location where to print debugging values
162   1020 FORMAT(A,1P4E11.3)  #include "THSICE_DEBUG.h"
163    
164        dt  = thSIce_deltaT   1010 FORMAT(A,I3,3F11.6)
165        dBug = .FALSE.   1020 FORMAT(A,1P4E14.6)
166  c     dBug = ( bi.EQ.3 .AND. i.EQ.15 .AND. j.EQ.11 )  
167  c     dBug = ( bi.EQ.6 .AND. i.EQ.18 .AND. j.EQ.10 )  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
168        IF (dBug) WRITE(6,'(A,2I4,2I2)') 'ThSI_SOLVE4T: i,j=',i,j,bi,bj  
169    #ifdef ALLOW_AUTODIFF_TAMC
170        if ( hi.lt.himin ) then        act1 = bi - myBxLo(myThid)
171  C If hi < himin, melt the ice.        max1 = myBxHi(myThid) - myBxLo(myThid) + 1
172           STOP 'THSICE_SOLVE4TEMP: should not enter if hi<himin'        act2 = bj - myByLo(myThid)
173        endif        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  devdt        = comlev1_thsice_1, key=ikey_1
195    CADJ STORE  df0dt        = comlev1_thsice_1, key=ikey_1
196    CADJ STORE  flxexceptsw  = comlev1_thsice_1, key=ikey_1
197    CADJ STORE  flxsw(i,j)   = comlev1_thsice_1, key=ikey_1
198    CADJ STORE  qic1(i,j)    = comlev1_thsice_1, key=ikey_1
199    CADJ STORE  qic2(i,j)    = comlev1_thsice_1, key=ikey_1
200    CADJ STORE  tsrf(i,j)    = comlev1_thsice_1, key=ikey_1
201    #endif
202            IF ( iceMask(i,j).GT.0. _d 0) THEN
203              hi      = hIce(i,j)
204              hs      = hSnow(i,j)
205              Tf      = tFrz(i,j)
206              netSW   = flxSW(i,j)
207              Tsf     = tSrf(i,j)
208              qicen(1)= qIc1(i,j)
209              qicen(2)= qIc2(i,j)
210    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
211    #ifdef ALLOW_DBUG_THSICE
212              IF ( dBug(i,j,bi,bj) ) WRITE(6,'(A,2I4,2I2)')
213         &         'ThSI_SOLVE4T: i,j=',i,j,bi,bj
214    #endif
215              IF ( hi.LT.hIceMin ) THEN
216    C If hi < hIceMin, melt the ice.
217                 STOP 'THSICE_SOLVE4TEMP: should not enter if hi<hIceMin'
218              ENDIF
219    
220  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
221    
222  C fractional snow cover  C fractional snow cover
223        frsnow = 0. _d 0        frsnow = 0. _d 0
224        if (hs .gt. 0. _d 0) frsnow = 1. _d 0        IF (hs .GT. 0. _d 0) frsnow = 1. _d 0
225    
226  C Compute SW flux absorbed at surface and penetrating to layer 1.  C Compute SW flux absorbed at surface and penetrating to layer 1.
227        fswpen  = flxSW * (1. _d 0 - frsnow) * i0        fswpen  = netSW * (1. _d 0 - frsnow) * i0swFrac
228        fswocn = fswpen * exp(-ksolar*hi)        fswocn = fswpen * exp(-ksolar*hi)
229        fswint = fswpen - fswocn        fswint = fswpen - fswocn
230    
231        fswdn = flxSW - fswpen        fswdn = netSW - fswpen
232    
233  C Compute conductivity terms at layer interfaces.  C Compute conductivity terms at layer interfaces.
234    
235        k12 = 4. _d 0*kice*ksnow / (ksnow*hi + 4. _d 0*kice*hs)        k12 = 4. _d 0*kIce*kSnow / (kSnow*hi + 4. _d 0*kIce*hs)
236        k32 = 2. _d 0*kice  / hi        k32 = 2. _d 0*kIce  / hi
237    
238  C compute ice temperatures  C compute ice temperatures
239        a1 = cpice        a1 = cpIce
240        b1 = qicen(1) + (cpwater-cpice )*Tmlt1 - Lfresh        b1 = qicen(1) + (cpWater-cpIce )*Tmlt1 - Lfresh
241        c1 = Lfresh * Tmlt1        c1 = Lfresh * Tmlt1
242        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
243        Tice(2) = (Lfresh-qicen(2)) / cpice        Tice(2) = (Lfresh-qicen(2)) / cpIce
244    
245        if (Tice(1).gt.0. _d 0 .or. Tice(2).gt.0. _d 0) then        IF (Tice(1).GT.0. _d 0 ) THEN
246            write (6,*) 'BBerr Tice(1) > 0 = ',Tice(1)         WRITE (standardMessageUnit,'(A,I12,1PE14.6)')
247            write (6,*) 'BBerr Tice(2) > 0 = ',Tice(2)       &       ' BBerr: Tice(1) > 0 ; it=', myIter, qicen(1)
248        endif         WRITE (standardMessageUnit,'(A,4I5,2F11.4)')
249        IF (dBug) WRITE(6,1010) 'ThSI_SOLVE4T: k, Ts, Tice=',0,Tsf,Tice       &      ' BBerr: i,j,bi,bj,Tice = ',i,j,bi,bj,Tice
250          ENDIF
251          IF ( Tice(2).GT.0. _d 0) THEN
252           WRITE (standardMessageUnit,'(A,I12,1PE14.6)')
253         &       ' BBerr: Tice(2) > 0 ; it=', myIter, qicen(2)
254           WRITE (standardMessageUnit,'(A,4I5,2F11.4)')
255         &      ' BBerr: i,j,bi,bj,Tice = ',i,j,bi,bj,Tice
256          ENDIF
257    #ifdef ALLOW_DBUG_THSICE
258          IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)
259         &  'ThSI_SOLVE4T: k, Ts, Tice=',0,Tsf,Tice
260    #endif
261    
262  C Compute coefficients used in quadratic formula.  C Compute coefficients used in quadratic formula.
263    
264        a10 = rhoi*cpice *hi/(2. _d 0*dt) +        a10 = rhoi*cpIce *hi/(2. _d 0*dt) +
265       &      k32 * (4. _d 0*dt*k32 + rhoi*cpice *hi)       &      k32 * (4. _d 0*dt*k32 + rhoi*cpIce *hi)
266       &       / (6. _d 0*dt*k32 + rhoi*cpice *hi)       &       / (6. _d 0*dt*k32 + rhoi*cpIce *hi)
267        b10 = -hi*        b10 = -hi*
268       &      (rhoi*cpice*Tice(1)+rhoi*Lfresh*Tmlt1/Tice(1))       &      (rhoi*cpIce*Tice(1)+rhoi*Lfresh*Tmlt1/Tice(1))
269       &       /(2. _d 0*dt)       &       /(2. _d 0*dt)
270       &      - k32 * (4. _d 0*dt*k32*Tf+rhoi*cpice *hi*Tice(2))       &      - k32 * (4. _d 0*dt*k32*Tf+rhoi*cpIce *hi*Tice(2))
271       &       / (6. _d 0*dt*k32 + rhoi*cpice *hi) - fswint       &       / (6. _d 0*dt*k32 + rhoi*cpIce *hi) - fswint
272        c1 = rhoi*Lfresh*hi*Tmlt1 / (2. _d 0*dt)        c1 = rhoi*Lfresh*hi*Tmlt1 / (2. _d 0*dt)
273    
274  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
275  C Compute new surface and internal temperatures; iterate until  C Compute new surface and internal temperatures; iterate until
276  C Tsfc converges.  C Tsfc converges.
277    
278          IF ( useBlkFlx ) THEN
279            iterMax = nitMaxTsf
280          ELSE
281            iterMax = 1
282          ENDIF
283          dTsf = Terrmax
284    
285  C ----- begin iteration  -----  C ----- begin iteration  -----
286        do 100 k = 1,nitMaxTsf        DO k = 1,iterMax
287    
288    #ifdef ALLOW_AUTODIFF_TAMC
289           ikey_3 = (ikey_1-1)*MaxTsf + k
290    #endif
291    
292    #ifdef ALLOW_AUTODIFF_TAMC
293    CADJ STORE tsf          = comlev1_thsice_3, key=ikey_3
294    CADJ STORE dtsf         = comlev1_thsice_3, key=ikey_3
295    CADJ STORE df0dt        = comlev1_thsice_3, key=ikey_3
296    CADJ STORE flxexceptsw  = comlev1_thsice_3, key=ikey_3
297    #endif
298           IF ( ABS(dTsf).GE.Terrmax ) THEN
299    
300  C Save temperatures at start of iteration.  C Save temperatures at start of iteration.
301  c        Tsf_start = Tsf  c        Tsf_start = Tsf
302    
303          IF ( useBlkFlx ) THEN          IF ( useBlkFlx ) THEN
304    #ifdef ALLOW_AUTODIFF_TAMC
305    CADJ STORE tsf          = comlev1_thsice_3, key=ikey_3
306    #endif
307  C Compute top surface flux.  C Compute top surface flux.
308           if (hs.gt.3. _d -1) then           IF (hs.GT.3. _d -1) THEN
309                iceornot=2                iceornot=2
310           else           ELSE
311                iceornot=1                iceornot=1
312           endif           ENDIF
313           CALL THSICE_GET_BULKF(           IF ( useBulkForce ) THEN
314       I                        iceornot, Tsf,            CALL THSICE_GET_BULKF(
315       O                        flxExceptSw, df0dT, evap, dEvdT,       I                         iceornot, Tsf,
316       I                        i,j,bi,bj,myThid )       O                         flxExceptSw, df0dT, evap, dEvdT,
317           flx0 = fswdn + flxExceptSw       I                         i,j,bi,bj,myThid )
318             ELSEIF ( useEXF ) THEN
319              CALL THSICE_GET_EXF  (
320         I                         iceornot, Tsf,
321         O                         flxExceptSw, df0dT, evap, dEvdT,
322         I                         i,j,bi,bj,myThid )
323             ENDIF
324          ELSE          ELSE
325           flx0 = fswdn + flxExcSw(1)           flxExceptSw = flxExSW(i,j,1)
326           df0dT = flxExcSw(2)           df0dT = flxExSW(i,j,2)
327          ENDIF          ENDIF
328          IF ( dBug ) WRITE(6,1020) 'ThSI_SOLVE4T: flx0,df0dT,k12,D=',          flx0 = fswdn + flxExceptSw
329       &                                     flx0,df0dT,k12,k12-df0dT  #ifdef ALLOW_DBUG_THSICE
330          IF ( dBug(i,j,bi,bj) ) WRITE(6,1020)
331         &  'ThSI_SOLVE4T: flx0,df0dT,k12,D=', flx0,df0dT,k12,k12-df0dT
332    #endif
333    
334  C Compute new top layer and surface temperatures.  C Compute new top layer and surface temperatures.
335  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
336  C with different coefficients.  C with different coefficients.
337    
338    #ifdef ALLOW_AUTODIFF_TAMC
339    CADJ STORE tsf  = comlev1_thsice_3, key=ikey_3
340    #endif
341           a1 = a10 - k12*df0dT / (k12-df0dT)           a1 = a10 - k12*df0dT / (k12-df0dT)
342           b1 = b10 - k12*(flx0-df0dT*Tsf) / (k12-df0dT)           b1 = b10 - k12*(flx0-df0dT*Tsf) / (k12-df0dT)
343           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)
344           dTsf = (flx0 + k12*(Tice(1)-Tsf)) / (k12-df0dT)           dTsf = (flx0 + k12*(Tice(1)-Tsf)) / (k12-df0dT)
345           Tsf = Tsf + dTsf           Tsf = Tsf + dTsf
346           if (Tsf .gt. 0. _d 0) then           IF (Tsf .GT. 0. _d 0) THEN
347              IF(dBug) WRITE(6,1010) 'ThSI_SOLVE4T: k,ts,t1,dTs=',  #ifdef ALLOW_DBUG_THSICE
348       &                                          k,Tsf,Tice(1),dTsf              IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)
349         &        'ThSI_SOLVE4T: k,ts,t1,dTs=', k,Tsf,Tice(1),dTsf
350    #endif
351              a1 = a10 + k12              a1 = a10 + k12
352              b1 = b10          ! note b1 = b10 - k12*Tf0              b1 = b10          ! note b1 = b10 - k12*Tf0
353              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)
354              Tsf = 0. _d 0              Tsf = 0. _d 0
355             IF ( useBlkFlx ) THEN             IF ( useBlkFlx ) THEN
356              if (hs.gt.3. _d -1) then              IF (hs.GT.3. _d -1) THEN
357                   iceornot=2                   iceornot=2
358              else              ELSE
359                   iceornot=1                   iceornot=1
360              endif              ENDIF
361              CALL THSICE_GET_BULKF(              IF ( useBulkForce ) THEN
362       I                        iceornot, Tsf,               CALL THSICE_GET_BULKF(
363       O                        flxExceptSw, df0dT, evap, dEvdT,       I                            iceornot, Tsf,
364       I                        i,j,bi,bj,myThid )       O                            flxExceptSw, df0dT, evap, dEvdT,
365              flx0 = fswdn + flxExceptSw       I                            i,j,bi,bj,myThid )
366                ELSEIF ( useEXF ) THEN
367                 CALL THSICE_GET_EXF  (
368         I                            iceornot, Tsf,
369         O                            flxExceptSw, df0dT, evap, dEvdT,
370         I                            i,j,bi,bj,myThid )
371                ENDIF
372              dTsf = 0. _d 0              dTsf = 0. _d 0
373             ELSE             ELSE
374              flx0 = fswdn + flxExcSw(0)              flxExceptSw = flxExSW(i,j,0)
375              dTsf = 1000.              dTsf = 1000.
376              df0dT = 0.              df0dT = 0.
377             ENDIF             ENDIF
378              Tsf  = 0. _d 0             flx0 = fswdn + flxExceptSw
379           endif           ENDIF
380    
381  C Check for convergence.  If no convergence, then repeat.  C Check for convergence.  If no convergence, then repeat.
382  C  C
383  C Convergence test: Make sure Tsfc has converged, within prescribed error.    C Convergence test: Make sure Tsfc has converged, within prescribed error.
384  C (Energy conservation is guaranteed within machine roundoff, even  C (Energy conservation is guaranteed within machine roundoff, even
385  C if Tsfc has not converged.)  C if Tsfc has not converged.)
386  C If no convergence, then repeat.  C If no convergence, then repeat.
387    
388           IF ( dBug ) WRITE(6,1010) 'ThSI_SOLVE4T: k,ts,t1,dTs=',  #ifdef ALLOW_DBUG_THSICE
389       &                                            k,Tsf,Tice(1),dTsf           IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)
390           IF ( useBlkFlx ) THEN       &     'ThSI_SOLVE4T: k,ts,t1,dTs=', k,Tsf,Tice(1),dTsf
391            if (abs(dTsf).lt.Terrmax) then  #endif
392              goto 150           IF ( useBlkFlx .AND. k.EQ.nitMaxTsf
393            elseif (k.eq.nitMaxTsf) then       &                  .AND. ABS(dTsf).GE.Terrmax ) THEN
394              write (6,*) 'BB: thermw conv err ',i,j,bi,bj,dTsf              WRITE (6,'(A,4I4,I12,F15.9)')
395              write (6,*) 'BB: thermw conv err, iceheight ', hi       &                 ' BB: not converge: i,j,it,hi=',i,j,bi,bj,
396              write (6,*) 'BB: thermw conv err: Tsf, flx0', Tsf,flx0       &                                                    myIter,hi
397              if (Tsf.lt.-70. _d 0) stop  !QQQQ fails              WRITE (6,*) 'BB: not converge: Tsf, dTsf=', Tsf,dTsf
398            endif              WRITE (6,*) 'BB: not converge: flx0,dfdT=',flx0,df0dT
399           ELSE              IF (Tsf.LT.-70. _d 0) STOP
             goto 150  
400           ENDIF           ENDIF
401    
402  100   continue  ! surface temperature iteration         ENDIF
403  150   continue        ENDDO
404  C ------ end iteration ------------  C ------ end iteration ------------
405    
 c        if (Tsf.lt.-70. _d 0) then  
 cQQ        print*,'QQQQQ Tsf =',Tsf  
 cQQ        stop  !QQQQ fails  
 c        endif  
   
406  C Compute new bottom layer temperature.  C Compute new bottom layer temperature.
407    
408    #ifdef ALLOW_AUTODIFF_TAMC
409    CADJ STORE  Tice(:)      = comlev1_thsice_1, key=ikey_1
410    CADJ STORE  df0dt        = comlev1_thsice_1, key=ikey_1
411    #endif
412        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)
413       &        + rhoi*cpice *hi*Tice(2))       &        + rhoi*cpIce *hi*Tice(2))
414       &         /(6. _d 0*dt*k32 + rhoi*cpice *hi)       &         /(6. _d 0*dt*k32 + rhoi*cpIce *hi)
415        IF (dBug) WRITE(6,1010) 'ThSI_SOLVE4T: k, Ts, Tice=',k,Tsf,Tice  #ifdef ALLOW_DBUG_THSICE
416          IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)
417         &  'ThSI_SOLVE4T: k, Ts, Tice=',k,Tsf,Tice
418    #endif
419    
420  C Compute final flux values at surfaces.  C Compute final flux values at surfaces.
421    
422        fct    = k12*(Tsf-Tice(1))        fct    = k12*(Tsf-Tice(1))
423        flxCnB = 4. _d 0*kice *(Tice(2)-Tf)/hi        flxCnB(i,j) = 4. _d 0*kIce *(Tice(2)-Tf)/hi
424        flx0   = flx0 + df0dT*dTsf        flx0   = flx0 + df0dT*dTsf
425        IF ( useBlkFlx ) THEN        IF ( useBlkFlx ) THEN
         evpAtm = evap  
426  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
427          evpAtm = evap + dEvdT*dTsf          evpAtm(i,j) = 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  
428        ELSE        ELSE
429          flxAtm = 0.  C- WARNING: Evap & +Evap*Lfresh are missing ! (but only affects Diagnostics)
430          evpAtm = 0.          evpAtm(i,j) = 0.
431        ENDIF        ENDIF
432        sHeating = flx0 - fct  C-    energy flux to Atmos: use net short-wave flux at surf. and
433    C     use latent heat = Lvap (energy=0 for liq. water at 0.oC)
434          flxAtm(i,j) = netSW + flxExceptSw
435         &                    + df0dT*dTsf + evpAtm(i,j)*Lfresh
436    C-    excess of energy @ surface (used for surface melting):
437          sHeat(i,j) = flx0 - fct
438    
439  C-    SW flux at sea-ice base left to the ocean  C-    SW flux at sea-ice base left to the ocean
440        flxSW = fswocn        flxSW(i,j) = fswocn
441    
442        IF (dBug) WRITE(6,1020) 'ThSI_SOLVE4T: flx0,fct,Dif,flxCnB=',  #ifdef ALLOW_DBUG_THSICE
443       &    flx0,fct,flx0-fct,flxCnB        IF ( dBug(i,j,bi,bj) ) WRITE(6,1020)
444         &  'ThSI_SOLVE4T: flx0,fct,Dif,flxCnB=',
445         &                 flx0,fct,flx0-fct,flxCnB(i,j)
446    #endif
447    
448  C Compute new enthalpy for each layer.  C Compute new enthalpy for each layer.
449    
450        qicen(1) = -cpwater*Tmlt1 + cpice *(Tmlt1-Tice(1)) +        qicen(1) = -cpWater*Tmlt1 + cpIce *(Tmlt1-Tice(1))
451       &              Lfresh*(1. _d 0-Tmlt1/Tice(1))       &            + Lfresh*(1. _d 0-Tmlt1/Tice(1))
452        qicen(2) = -cpice *Tice(2) + Lfresh        qicen(2) = -cpIce *Tice(2) + Lfresh
453    
454  C Make sure internal ice temperatures do not exceed Tmlt.  C Make sure internal ice temperatures do not exceed Tmlt.
455  C (This should not happen for reasonable values of i0.)  C (This should not happen for reasonable values of i0swFrac)
456    
457        if (Tice(1) .ge. Tmlt1) then        IF (Tice(1) .GE. Tmlt1) THEN
458          write (6,'(A,2I4,2I3,1P2E14.6)')          WRITE (6,'(A,2I4,2I3,1P2E14.6)')
459       &   'BBerr - Bug: IceT(1) > Tmlt',i,j,bi,bj,Tice(1),Tmlt1       &   ' BBerr - Bug: IceT(1) > Tmlt',i,j,bi,bj,Tice(1),Tmlt1
460        endif        ENDIF
461        if (Tice(2) .ge. 0. _d 0) then        IF (Tice(2) .GE. 0. _d 0) THEN
462         write (6,'(A,2I4,2I3,1P2E14.6)')         WRITE (6,'(A,2I4,2I3,1P2E14.6)')
463       &   'BBerr - Bug: IceT(2) > 0',i,j,bi,bj,Tice(2)       &   ' BBerr - Bug: IceT(2) > 0',i,j,bi,bj,Tice(2)
464        endif        ENDIF
465    
466    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
467    C--    Update Sea-Ice state :
468              tSrf(i,j)  = Tsf
469              tIc1(i,j)  = Tice(1)
470              tic2(i,j)  = Tice(2)
471              qIc1(i,j)  = qicen(1)
472              qIc2(i,j)  = qicen(2)
473    c         dTsrf(i,j) = dTsf
474              IF ( .NOT.useBlkFlx ) dTsrf(i,j) = dTsf
475    c         sHeat(i,j) = sHeating
476    c         flxCnB(i,j)= flxCnB
477    c         flxAtm(i,j)= flxAtm
478    c         evpAtm(i,j)= evpAtm
479    #ifdef ALLOW_DBUG_THSICE
480              IF ( dBug(i,j,bi,bj) ) THEN
481               WRITE(6,1020) 'ThSI_SOLV_4T: Tsf, Tice(1,2), dTsurf=',
482         &                              Tsf, Tice, dTsf
483               WRITE(6,1020) 'ThSI_SOLV_4T: sHeat, flxCndBt, Qice =',
484         &           sHeat(i,j), flxCnB(i,j), qicen
485               WRITE(6,1020) 'ThSI_SOLV_4T: flxA, evpA, fxSW_bf,af=',
486         &           flxAtm(i,j), evpAtm(i,j), netSW, flxSW(i,j)
487              ENDIF
488    #endif
489            ELSE
490              IF ( .NOT.useBlkFlx ) dTsrf(i,j) = 0. _d 0
491            ENDIF
492           ENDDO
493          ENDDO
494  #endif  /* ALLOW_THSICE */  #endif  /* ALLOW_THSICE */
495    
496  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.16
  ViewVC Help
Powered by ViewVC 1.1.22