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