/[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.17 by jmc, Tue May 1 02:26:50 2007 UTC revision 1.22 by mlosch, Fri Oct 15 15:15:31 2010 UTC
# Line 27  C ADAPTED FROM: Line 27  C ADAPTED FROM:
27  C LANL CICE.v2.0.2  C LANL CICE.v2.0.2
28  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
29  C.. thermodynamics (vertical physics) based on M. Winton 3-layer model  C.. thermodynamics (vertical physics) based on M. Winton 3-layer model
30  C.. See Bitz, C. M. and W. H. Lipscomb, 1999:  "An energy-conserving  C.. See Bitz, C. M. and W. H. Lipscomb, 1999:  An energy-conserving
31  C..       thermodynamic sea ice model for climate study."  J. Geophys.  C..       thermodynamic sea ice model for climate study.
32  C..       Res., 104, 15669 - 15677.  C..       J. Geophys. Res., 104, 15669 - 15677.
33  C..     Winton, M., 1999:  "A reformulated three-layer sea ice model."  C..     Winton, M., 1999:  "A reformulated three-layer sea ice model."
34  C..       Submitted to J. Atmos. Ocean. Technol.  C..       Submitted to J. Atmos. Ocean. Technol.
35  C.. authors Elizabeth C. Hunke and William Lipscomb  C.. authors Elizabeth C. Hunke and William Lipscomb
# Line 44  C     !USES: Line 44  C     !USES:
44    
45  C     == Global variables ===  C     == Global variables ===
46  #include "EEPARAMS.h"  #include "EEPARAMS.h"
47    #include "SIZE.h"
48  #include "THSICE_SIZE.h"  #include "THSICE_SIZE.h"
49  #include "THSICE_PARAMS.h"  #include "THSICE_PARAMS.h"
50  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
 # include "SIZE.h"  
51  # include "tamc.h"  # include "tamc.h"
52  # include "tamc_keys.h"  # include "tamc_keys.h"
53  #endif  #endif
# Line 63  C     dBugFlag    :: allow to print debu Line 63  C     dBugFlag    :: allow to print debu
63  C     useBulkForce:: use surf. fluxes from bulk-forcing external S/R  C     useBulkForce:: use surf. fluxes from bulk-forcing external S/R
64  C     useEXF      :: use surf. fluxes from exf          external S/R  C     useEXF      :: use surf. fluxes from exf          external S/R
65  C---  Input:  C---  Input:
66  C         iceMask :: sea-ice fractional mask [0-1]  C     iceMask     :: sea-ice fractional mask [0-1]
67  C  hIce    (hi)   :: ice height [m]  C     hIce        :: ice height [m]
68  C  hSnow   (hs)   :: snow height [m]  C     hSnow       :: snow height [m]
69  C  tFrz    (Tf)   :: sea-water freezing temperature [oC] (function of S)  C     tFrz        :: sea-water freezing temperature [oC] (function of S)
70  C  flxExSW  (=)   :: surf. heat flux (+=down) except SW, function of surf. temp Ts:  C     flxExSW     :: surf. heat flux (+=down) except SW, function of surf. temp Ts:
71  C                    0: Flx(Ts=0) ; 1: Flx(Ts=Ts^n) ; 2: d.Flx/dTs(Ts=Ts^n)  C                    0: Flx(Ts=0) ; 1: Flx(Ts=Ts^n) ; 2: d.Flx/dTs(Ts=Ts^n)
72  C---  Modified (input&output):  C---  Modified (input&output):
73  C  flxSW  (netSW) :: net Short-Wave flux (+=down) [W/m2]: input= at surface  C     flxSW       :: net Short-Wave flux (+=down) [W/m2]: input= at surface
74  C           (=)   ::                 output= below sea-ice, into the ocean  C                 ::                 output= below sea-ice, into the ocean
75  C  tSrf    (Tsf)  :: surface (ice or snow) temperature  C     tSrf        :: surface (ice or snow) temperature
76  C  qIc1   (qicen) :: ice enthalpy (J/kg), 1rst level  C     qIc1        :: ice enthalpy (J/kg), 1rst level
77  C  qIc2   (qicen) :: ice enthalpy (J/kg), 2nd level  C     qIc2        :: ice enthalpy (J/kg), 2nd level
78  C---  Output  C---  Output
79  C  tIc1    (Tice) :: temperature of ice layer 1 [oC]  C     tIc1        :: temperature of ice layer 1 [oC]
80  C  tIc2    (Tice) :: temperature of ice layer 2 [oC]  C     tIc2        :: temperature of ice layer 2 [oC]
81  C  dTsrf   (dTsf) :: surf. temp adjusment: Ts^n+1 - Ts^n  C     dTsrf       :: surf. temp adjusment: Ts^n+1 - Ts^n
82  C  sHeat(sHeating):: surf heating flux left to melt snow or ice (= Atmos-conduction)  C     sHeat       :: surf heating flux left to melt snow or ice (= Atmos-conduction)
83  C  flxCnB   (=)   :: heat flux conducted through the ice to bottom surface  C     flxCnB      :: heat flux conducted through the ice to bottom surface
84  C  flxAtm   (=)   :: net flux of energy from the atmosphere [W/m2] (+=down)  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)  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)  C     evpAtm      :: evaporation to the atmosphere [kg/m2/s] (>0 if evaporate)
87  C---  Input:  C---  Input:
88  C     myTime      :: current Time of simulation [s]  C     myTime      :: current Time of simulation [s]
89  C     myIter      :: current Iteration number in simulation  C     myIter      :: current Iteration number in simulation
# Line 106  C     myThid      :: my Thread Id number Line 106  C     myThid      :: my Thread Id number
106        _RL qIc2   (siLo:siHi,sjLo:sjHi)        _RL qIc2   (siLo:siHi,sjLo:sjHi)
107        _RL tIc1   (siLo:siHi,sjLo:sjHi)        _RL tIc1   (siLo:siHi,sjLo:sjHi)
108        _RL tIc2   (siLo:siHi,sjLo:sjHi)        _RL tIc2   (siLo:siHi,sjLo:sjHi)
 c     _RL dTsrf  (siLo:siHi,sjLo:sjHi)  
109        _RL dTsrf  (iMin:iMax,jMin:jMax)        _RL dTsrf  (iMin:iMax,jMin:jMax)
110        _RL sHeat  (siLo:siHi,sjLo:sjHi)        _RL sHeat  (siLo:siHi,sjLo:sjHi)
111        _RL flxCnB (siLo:siHi,sjLo:sjHi)        _RL flxCnB (siLo:siHi,sjLo:sjHi)
# Line 119  CEOP Line 118  CEOP
118    
119  #ifdef ALLOW_THSICE  #ifdef ALLOW_THSICE
120  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
 C---  local copy of input/output argument list variables (see description above)  
 c     _RL flxExcSw(0:2)  
       _RL Tf  
       _RL hi  
       _RL hs  
       _RL netSW  
       _RL Tsf  
       _RL qicen(nlyr)  
       _RL Tice (nlyr)  
 c     _RL sHeating  
 c     _RL flxCnB  
       _RL dTsf  
 c     _RL flxAtm  
 c     _RL evpAtm  
   
121  C     == Local Variables ==  C     == Local Variables ==
122    C     useBlkFlx    :: use some bulk-formulae to compute fluxes over sea-ice
123    C                  :: (otherwise, fluxes are passed as argument from AIM)
124    C     iterate4Tsf  :: to stop to iterate when all icy grid pts Tsf did converged
125    C     iceFlag      :: True= do iterate for Surf.Temp ; False= do nothing
126    C     frsnow       :: fractional snow cover
127    C     fswpen       :: SW penetrating beneath surface (W m-2)
128    C     fswdn        :: SW absorbed at surface (W m-2)
129    C     fswint       :: SW absorbed in ice (W m-2)
130    C     fswocn       :: SW passed through ice to ocean (W m-2)
131    C     Tsf          :: surface (ice or snow) temperature   (local copy of tSrf)
132    C     flx0exSW     :: net surface heat flux over melting snow/ice, except short-wave.
133    C     flxTexSW     :: net surface heat flux, except short-wave (W/m2)
134    C     dFlxdT       :: deriv of flxNet wrt Tsf (W m-2 deg-1)
135    C     evap0        :: evaporation over melting snow/ice [kg/m2/s] (>0 if evaporate)
136    C     evapT        :: evaporation over snow/ice [kg/m2/s] (>0 if evaporate)
137    C     dEvdT        :: derivative of evap. with respect to Tsf [kg/m2/s/K]
138    C     flxNet       :: net surf heat flux (+=down), from Atmos. to sea-ice (W m-2)
139    C     netSW        :: net Short-Wave flux at surface (+=down) [W/m2]
140    C     fct          :: heat conducted to top surface
141    C     k12, k32     :: thermal conductivity terms
142    C     a10,b10,c10  :: coefficients in quadratic eqn for T1
143    C     a1, b1, c1   :: coefficients in quadratic eqn for T1
144    C     dt           :: timestep
145          LOGICAL useBlkFlx
146          LOGICAL iterate4Tsf
147          LOGICAL iceFlag(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
148        INTEGER i,j        INTEGER i,j
149        INTEGER k, iterMax        INTEGER k, iterMax, ii, jj, icount
150          _RL  frsnow
151        _RL  frsnow        ! fractional snow cover        _RL  fswpen
152        _RL  fswpen        ! SW penetrating beneath surface (W m-2)        _RL  fswdn
153        _RL  fswdn         ! SW absorbed at surface (W m-2)        _RL  fswint
154        _RL  fswint        ! SW absorbed in ice (W m-2)        _RL  fswocn
155        _RL  fswocn        ! SW passed through ice to ocean (W m-2)        _RL  Tsf     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
156        _RL  flxExceptSw   ! net surface heat flux, except short-wave (W/m2)        _RL  flx0exSW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
157  C     evap           :: evaporation over snow/ice [kg/m2/s] (>0 if evaporate)        _RL  flxTexSW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
158  C     dEvdT          :: derivative of evap. with respect to Tsf [kg/m2/s/K]        _RL  dFlxdT  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
159        _RL  evap, dEvdT        _RL  evap0   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
160        _RL  flx0          ! net surf heat flux, from Atmos. to sea-ice (W m-2)        _RL  evapT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
161        _RL  fct           ! heat conducted to top surface        _RL  dEvdT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
162        _RL  df0dT         ! deriv of flx0 wrt Tsf (W m-2 deg-1)        _RL  flxNet
163        _RL  k12, k32      ! thermal conductivity terms        _RL  fct
164        _RL  a10, b10      ! coefficients in quadratic eqn for T1        _RL  k12     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
165        _RL  a1, b1, c1    ! coefficients in quadratic eqn for T1        _RL  k32
166  c     _RL  Tsf_start     ! old value of Tsf        _RL  a10     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
167        _RL  dt            ! timestep        _RL  b10     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
168          _RL  c10     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
169          _RL  a1, b1, c1
170          _RL  dt
171        _RL  recip_dhSnowLin        _RL  recip_dhSnowLin
172        INTEGER iceornot  #ifdef ALLOW_DBUG_THSICE
173        LOGICAL useBlkFlx        _RL  netSW
174    #endif
175    
176  C-    define grid-point location where to print debugging values  C-    Define grid-point location where to print debugging values
177  #include "THSICE_DEBUG.h"  #include "THSICE_DEBUG.h"
178    
179   1010 FORMAT(A,I3,3F11.6)   1010 FORMAT(A,I3,3F11.6)
# Line 168  C-    define grid-point location where t Line 182  C-    define grid-point location where t
182  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
183    
184  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
185        act1 = bi - myBxLo(myThid)  c     act1 = bi - myBxLo(myThid)
186        max1 = myBxHi(myThid) - myBxLo(myThid) + 1  c     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
187        act2 = bj - myByLo(myThid)  c     act2 = bj - myByLo(myThid)
188        max2 = myByHi(myThid) - myByLo(myThid) + 1  c     max2 = myByHi(myThid) - myByLo(myThid) + 1
189        act3 = myThid - 1  c     act3 = myThid - 1
190        max3 = nTx*nTy  c     max3 = nTx*nTy
191        act4 = ikey_dynamics - 1  c     act4 = ikey_dynamics - 1
192  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
193    
194        useBlkFlx = useEXF .OR. useBulkForce        useBlkFlx = useEXF .OR. useBulkForce
195          dt  = thSIce_dtTemp
196        IF ( dhSnowLin.GT.0. _d 0 ) THEN        IF ( dhSnowLin.GT.0. _d 0 ) THEN
197          recip_dhSnowLin = 1. _d 0 / dhSnowLin          recip_dhSnowLin = 1. _d 0 / dhSnowLin
198        ELSE        ELSE
199          recip_dhSnowLin = 0. _d 0          recip_dhSnowLin = 0. _d 0
200        ENDIF        ENDIF
201    
202        dt  = thSIce_dtTemp        iterate4Tsf = .FALSE.
203          icount = 0
204    C    
205        DO j = jMin, jMax        DO j = jMin, jMax
206         DO i = iMin, iMax         DO i = iMin, iMax
207  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
208            ikey_1 = i  c         ikey_1 = i
209       &         + sNx*(j-1)  c    &         + sNx*(j-1)
210       &         + sNx*sNy*act1  c    &         + sNx*sNy*act1
211       &         + sNx*sNy*max1*act2  c    &         + sNx*sNy*max1*act2
212       &         + sNx*sNy*max1*max2*act3  c    &         + sNx*sNy*max1*max2*act3
213       &         + sNx*sNy*max1*max2*max3*act4  c    &         + sNx*sNy*max1*max2*max3*act4
214  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
 C--  
215  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
216  CADJ STORE  devdt        = comlev1_thsice_1, key=ikey_1  cCADJ STORE  devdt        = comlev1_thsice_1, key=ikey_1
217  CADJ STORE  df0dt        = comlev1_thsice_1, key=ikey_1  cCADJ STORE  dFlxdT        = comlev1_thsice_1, key=ikey_1
218  CADJ STORE  flxexceptsw  = comlev1_thsice_1, key=ikey_1  cCADJ STORE  flxexceptsw  = comlev1_thsice_1, key=ikey_1
219  CADJ STORE  flxsw(i,j)   = comlev1_thsice_1, key=ikey_1  cCADJ STORE  flxsw(i,j)   = comlev1_thsice_1, key=ikey_1
220  CADJ STORE  qic1(i,j)    = comlev1_thsice_1, key=ikey_1  cCADJ STORE  qic1(i,j)    = comlev1_thsice_1, key=ikey_1
221  CADJ STORE  qic2(i,j)    = comlev1_thsice_1, key=ikey_1  cCADJ STORE  qic2(i,j)    = comlev1_thsice_1, key=ikey_1
222  CADJ STORE  tsrf(i,j)    = comlev1_thsice_1, key=ikey_1  cCADJ STORE  tsrf(i,j)    = comlev1_thsice_1, key=ikey_1
223  #endif  #endif
224          IF ( iceMask(i,j).GT.0. _d 0) THEN          IF ( iceMask(i,j).GT.0. _d 0) THEN
225            hi      = hIce(i,j)            iterate4Tsf  = .TRUE.
226            hs      = hSnow(i,j)            iceFlag(i,j) = .TRUE.
           Tf      = tFrz(i,j)  
           netSW   = flxSW(i,j)  
           Tsf     = tSrf(i,j)  
           qicen(1)= qIc1(i,j)  
           qicen(2)= qIc2(i,j)  
 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  
227  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
228            IF ( dBug(i,j,bi,bj) ) WRITE(6,'(A,2I4,2I2)')            IF ( dBug(i,j,bi,bj) ) WRITE(6,'(A,2I4,2I2)')
229       &         'ThSI_SOLVE4T: i,j=',i,j,bi,bj       &         'ThSI_SOLVE4T: i,j=',i,j,bi,bj
230  #endif  #endif
231            IF ( hi.LT.hIceMin ) THEN            IF ( hIce(i,j).LT.hIceMin ) THEN
232  C If hi < hIceMin, melt the ice.  C     if hi < hIceMin, melt the ice.
233               STOP 'THSICE_SOLVE4TEMP: should not enter if hi<hIceMin'  C     keep the position of this problem but do the stop later
234               ii = i
235               jj = j
236               icount = icount + 1
237            ENDIF            ENDIF
238    
239  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C--   Fractional snow cover:
240    C     assume a linear distribution of snow thickness, with dhSnowLin slope,
241    C      from hs-dhSnowLin to hs+dhSnowLin if full ice & snow cover.
242    C     frsnow = fraction of snow over the ice-covered part of the grid cell
243              IF ( hSnow(i,j) .GT. iceMask(i,j)*dhSnowLin ) THEN
244               frsnow = 1. _d 0
245              ELSE
246               frsnow = hSnow(i,j)*recip_dhSnowLin/iceMask(i,j)
247               IF ( frsnow.GT.0. _d 0 ) frsnow = SQRT(frsnow)
248              ENDIF
249    
250  C fractional snow cover  C--   Compute SW flux absorbed at surface and penetrating to layer 1.
251  C  assume a linear distribution of snow thickness, with dhSnowLin slope,            fswpen = flxSW(i,j) * (1. _d 0 - frsnow) * i0swFrac
252  C   from hs-dhSnowLin to hs+dhSnowLin if full ice & snow cover.            fswocn = fswpen * exp(-ksolar*hIce(i,j))
253  C frsnow = fraction of snow over the ice-covered part of the grid cell            fswint = fswpen - fswocn
254        IF ( hs .GT. iceMask(i,j)*dhSnowLin ) THEN            fswdn  = flxSW(i,j) - fswpen
255          frsnow = 1. _d 0  
256        ELSE  C     Initialise Atmospheric surf. heat flux with net SW flux at the surface
257          frsnow = hs*recip_dhSnowLin/iceMask(i,j)            flxAtm(i,j) = flxSW(i,j)
258          IF ( frsnow.GT.0. _d 0 ) frsnow = SQRT(frsnow)  C     SW flux at sea-ice base left to the ocean
259        ENDIF            flxSW(i,j) = fswocn
260    C     Initialise surface Heating with SW contribution
261              sHeat(i,j) = fswdn
262    
263    C--   Compute conductivity terms at layer interfaces.
264              k12(i,j) = 4. _d 0*kIce*kSnow
265         &             / (kSnow*hIce(i,j) + 4. _d 0*kIce*hSnow(i,j))
266              k32      = 2. _d 0*kIce  / hIce(i,j)
267    
268    C--   Compute ice temperatures
269              a1 = cpIce
270              b1 = qIc1(i,j) + (cpWater-cpIce )*Tmlt1 - Lfresh
271              c1 = Lfresh * Tmlt1
272              tIc1(i,j) = 0.5 _d 0 *(-b1 - SQRT(b1*b1-4. _d 0*a1*c1))/a1
273              tIc2(i,j) = (Lfresh-qIc2(i,j)) / cpIce
274    
 C Compute SW flux absorbed at surface and penetrating to layer 1.  
       fswpen  = netSW * (1. _d 0 - frsnow) * i0swFrac  
       fswocn = fswpen * exp(-ksolar*hi)  
       fswint = fswpen - fswocn  
   
       fswdn = netSW - fswpen  
   
 C Compute conductivity terms at layer interfaces.  
   
       k12 = 4. _d 0*kIce*kSnow / (kSnow*hi + 4. _d 0*kIce*hs)  
       k32 = 2. _d 0*kIce  / hi  
   
 C compute ice temperatures  
       a1 = cpIce  
       b1 = qicen(1) + (cpWater-cpIce )*Tmlt1 - Lfresh  
       c1 = Lfresh * Tmlt1  
       Tice(1) = 0.5 _d 0 *(-b1 - SQRT(b1*b1-4. _d 0*a1*c1))/a1  
       Tice(2) = (Lfresh-qicen(2)) / cpIce  
   
       IF (Tice(1).GT.0. _d 0 ) THEN  
        WRITE (standardMessageUnit,'(A,I12,1PE14.6)')  
      &       ' BBerr: Tice(1) > 0 ; it=', myIter, qicen(1)  
        WRITE (standardMessageUnit,'(A,4I5,2F11.4)')  
      &      ' BBerr: i,j,bi,bj,Tice = ',i,j,bi,bj,Tice  
       ENDIF  
       IF ( Tice(2).GT.0. _d 0) THEN  
        WRITE (standardMessageUnit,'(A,I12,1PE14.6)')  
      &       ' BBerr: Tice(2) > 0 ; it=', myIter, qicen(2)  
        WRITE (standardMessageUnit,'(A,4I5,2F11.4)')  
      &      ' BBerr: i,j,bi,bj,Tice = ',i,j,bi,bj,Tice  
       ENDIF  
275  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
276        IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)            IF (tIc1(i,j).GT.0. _d 0 ) THEN
277       &  'ThSI_SOLVE4T: k, Ts, Tice=',0,Tsf,Tice             WRITE (standardMessageUnit,'(A,I12,1PE14.6)')
278         &       ' BBerr: Tice(1) > 0 ; it=', myIter, qIc1(i,j)
279               WRITE (standardMessageUnit,'(A,4I5,2F11.4)')
280         &      ' BBerr: i,j,bi,bj,Tice = ',i,j,bi,bj,tIc1(i,j),tIc2(i,j)
281              ENDIF
282              IF ( tIc2(i,j).GT.0. _d 0) THEN
283               WRITE (standardMessageUnit,'(A,I12,1PE14.6)')
284         &       ' BBerr: Tice(2) > 0 ; it=', myIter, qIc2(i,j)
285               WRITE (standardMessageUnit,'(A,4I5,2F11.4)')
286         &      ' BBerr: i,j,bi,bj,Tice = ',i,j,bi,bj,tIc1(i,j),tIc2(i,j)
287              ENDIF
288              IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)
289         &      'ThSI_SOLVE4T: k, Ts, Tice=',0,tSrf(i,j),tIc1(i,j),tIc2(i,j)
290  #endif  #endif
291    
292  C Compute coefficients used in quadratic formula.  C--   Compute coefficients used in quadratic formula.
293    
294        a10 = rhoi*cpIce *hi/(2. _d 0*dt) +            a10(i,j) = rhoi*cpIce *hIce(i,j)/(2. _d 0*dt) +
295       &      k32 * (4. _d 0*dt*k32 + rhoi*cpIce *hi)       &          k32*( 4. _d 0*dt*k32 + rhoi*cpIce *hIce(i,j) )
296       &       / (6. _d 0*dt*k32 + rhoi*cpIce *hi)       &           / ( 6. _d 0*dt*k32 + rhoi*cpIce *hIce(i,j) )
297        b10 = -hi*            b10(i,j) = -hIce(i,j)*
298       &      (rhoi*cpIce*Tice(1)+rhoi*Lfresh*Tmlt1/Tice(1))       &          ( rhoi*cpIce*tIc1(i,j) + rhoi*Lfresh*Tmlt1/tIc1(i,j) )
299       &       /(2. _d 0*dt)       &           /(2. _d 0*dt)
300       &      - k32 * (4. _d 0*dt*k32*Tf+rhoi*cpIce *hi*Tice(2))       &        - k32*( 4. _d 0*dt*k32*tFrz(i,j)
301       &       / (6. _d 0*dt*k32 + rhoi*cpIce *hi) - fswint       &               +rhoi*cpIce*hIce(i,j)*tIc2(i,j) )
302        c1 = rhoi*Lfresh*hi*Tmlt1 / (2. _d 0*dt)       &           / ( 6. _d 0*dt*k32 + rhoi*cpIce *hIce(i,j) )
303         &        - fswint
304              c10(i,j) = rhoi*Lfresh*hIce(i,j)*Tmlt1 / (2. _d 0*dt)
305    
306            ELSE
307              iceFlag(i,j) = .FALSE.
308            ENDIF
309           ENDDO
310          ENDDO
311          IF ( icount .gt. 0 ) THEN
312           WRITE (standardMessageUnit,'(A,I5,A)')
313         &      'THSICE_SOLVE4TEMP: there are ',icount,
314         &      ' case(s) where hIce<hIceMin;'
315           WRITE (standardMessageUnit,'(A,I3,A1,I3,A)')
316         &      'THSICE_SOLVE4TEMP: the last one was at (',ii,',',jj,
317         &      ') with hIce = ', hIce(ii,jj)
318           STOP 'THSICE_SOLVE4TEMP: should not enter if hIce<hIceMin'
319          ENDIF
320    
321  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
322  C Compute new surface and internal temperatures; iterate until  C--   Get surface fluxes over melting surface
323  C Tsfc converges.  #ifdef ALLOW_AUTODIFF_TAMC
324          IF ( useBlkFlx ) THEN
325    #else
326          IF ( useBlkFlx .AND. iterate4Tsf  ) THEN
327    #endif
328            DO j = jMin, jMax
329             DO i = iMin, iMax
330               Tsf(i,j) = 0.
331             ENDDO
332            ENDDO
333            IF ( useEXF ) THEN
334               CALL THSICE_GET_EXF(   bi, bj, siLo, siHi, sjLo, sjHi,
335         I                            iMin, iMax, jMin, jMax,
336         I                            iceFlag, hSnow, Tsf,
337         O                            flx0exSW, dFlxdT, evap0, dEvdT,
338         I                            myTime, myIter, myThid )
339    C-    could add this "ifdef" to hide THSICE_GET_BULKF from TAF
340    c#ifdef ALLOW_BULK_FORCE
341            ELSEIF ( useBulkForce ) THEN
342               CALL THSICE_GET_BULKF( bi, bj, siLo, siHi, sjLo, sjHi,
343         I                            iMin, iMax, jMin, jMax,
344         I                            iceFlag, hSnow, Tsf,
345         O                            flx0exSW, dFlxdT, evap0, dEvdT,
346         I                            myTime, myIter, myThid )
347    c#endif /* ALLOW_BULK_FORCE */
348            ENDIF
349          ENDIF
350    
351    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
352    C---  Compute new surface and internal temperatures; iterate until
353    C     Tsfc converges.
354          DO j = jMin, jMax
355            DO i = iMin, iMax
356              Tsf(i,j)  = tSrf(i,j)
357              dTsrf(i,j) = Terrmax
358            ENDDO
359          ENDDO
360        IF ( useBlkFlx ) THEN        IF ( useBlkFlx ) THEN
361          iterMax = nitMaxTsf          iterMax = nitMaxTsf
362        ELSE        ELSE
363          iterMax = 1          iterMax = 1
364        ENDIF        ENDIF
       dTsf = Terrmax  
365    
366  C ----- begin iteration  -----  C ----- begin iteration  -----
367        DO k = 1,iterMax        DO k = 1,iterMax
368    #ifndef ALLOW_AUTODIFF_TAMC
369  #ifdef ALLOW_AUTODIFF_TAMC         IF ( iterate4Tsf ) THEN
370         ikey_3 = (ikey_1-1)*MaxTsf + k         iterate4Tsf = .FALSE.
371  #endif  #endif
372    
373  #ifdef ALLOW_AUTODIFF_TAMC         IF ( useBlkFlx ) THEN
374  CADJ STORE tsf          = comlev1_thsice_3, key=ikey_3  C--   Compute top surface flux.
375  CADJ STORE dtsf         = comlev1_thsice_3, key=ikey_3           IF ( useEXF ) THEN
376  CADJ STORE df0dt        = comlev1_thsice_3, key=ikey_3             CALL THSICE_GET_EXF(   bi, bj, siLo, siHi, sjLo, sjHi,
377  CADJ STORE flxexceptsw  = comlev1_thsice_3, key=ikey_3       I                            iMin, iMax, jMin, jMax,
378  #endif       I                            iceFlag, hSnow, Tsf,
379         IF ( ABS(dTsf).GE.Terrmax ) THEN       O                            flxTexSW, dFlxdT, evapT, dEvdT,
380         I                            myTime, myIter, myThid )
381  C Save temperatures at start of iteration.  C-    could add this "ifdef" to hide THSICE_GET_BULKF from TAF
382  c        Tsf_start = Tsf  c#ifdef ALLOW_BULK_FORCE
383             ELSEIF ( useBulkForce ) THEN
384          IF ( useBlkFlx ) THEN             CALL THSICE_GET_BULKF( bi, bj, siLo, siHi, sjLo, sjHi,
385  #ifdef ALLOW_AUTODIFF_TAMC       I                            iMin, iMax, jMin, jMax,
386  CADJ STORE tsf          = comlev1_thsice_3, key=ikey_3       I                            iceFlag, hSnow, Tsf,
387  #endif       O                            flxTexSW, dFlxdT, evapT, dEvdT,
388  C Compute top surface flux.       I                            myTime, myIter, myThid )
389           IF (hs.GT.3. _d -1) THEN  c#endif /* ALLOW_BULK_FORCE */
               iceornot=2  
          ELSE  
               iceornot=1  
          ENDIF  
          IF ( useBulkForce ) THEN  
           CALL THSICE_GET_BULKF(  
      I                         iceornot, Tsf,  
      O                         flxExceptSw, df0dT, evap, dEvdT,  
      I                         i,j,bi,bj,myThid )  
          ELSEIF ( useEXF ) THEN  
           CALL THSICE_GET_EXF  (  
      I                         iceornot, Tsf,  
      O                         flxExceptSw, df0dT, evap, dEvdT,  
      I                         i,j,bi,bj,myThid )  
390           ENDIF           ENDIF
391          ELSE         ELSE
392           flxExceptSw = flxExSW(i,j,1)           DO j = jMin, jMax
393           df0dT = flxExSW(i,j,2)            DO i = iMin, iMax
394          ENDIF             IF ( iceFlag(i,j) ) THEN
395          flx0 = fswdn + flxExceptSw               flxTexSW(i,j) = flxExSW(i,j,1)
396                 dFlxdT(i,j) = flxExSW(i,j,2)
397               ENDIF
398              ENDDO
399             ENDDO
400           ENDIF
401    
402    C--   Compute new top layer and surface temperatures.
403    C     If Tsfc is computed to be > 0 C, fix Tsfc = 0 and recompute T1
404    C     with different coefficients.
405           DO j = jMin, jMax
406            DO i = iMin, iMax
407             IF ( iceFlag(i,j) ) THEN
408              flxNet = sHeat(i,j) + flxTexSW(i,j)
409  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
410        IF ( dBug(i,j,bi,bj) ) WRITE(6,1020)            IF ( dBug(i,j,bi,bj) ) WRITE(6,1020)
411       &  'ThSI_SOLVE4T: flx0,df0dT,k12,D=', flx0,df0dT,k12,k12-df0dT       &     'ThSI_SOLVE4T: flxNet,dFlxdT,k12,D=',
412         &      flxNet, dFlxdT(i,j), k12(i,j), k12(i,j)-dFlxdT(i,j)
413  #endif  #endif
414    
 C Compute new top layer and surface temperatures.  
 C If Tsfc is computed to be > 0 C, fix Tsfc = 0 and recompute T1  
 C with different coefficients.  
   
415  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
416  CADJ STORE tsf  = comlev1_thsice_3, key=ikey_3  cCADJ STORE tsf  = comlev1_thsice_3, key=ikey_3
417  #endif  #endif
418           a1 = a10 - k12*df0dT / (k12-df0dT)            a1 = a10(i,j) - k12(i,j)*dFlxdT(i,j) / (k12(i,j)-dFlxdT(i,j))
419           b1 = b10 - k12*(flx0-df0dT*Tsf) / (k12-df0dT)            b1 = b10(i,j) - k12(i,j)*(flxNet-dFlxdT(i,j)*Tsf(i,j))
420           Tice(1) = -(b1 + SQRT(b1*b1-4. _d 0*a1*c1))/(2. _d 0*a1)       &                  /(k12(i,j)-dFlxdT(i,j))
421           dTsf = (flx0 + k12*(Tice(1)-Tsf)) / (k12-df0dT)            c1 = c10(i,j)
422           Tsf = Tsf + dTsf            tIc1(i,j)  = -(b1 + SQRT(b1*b1-4. _d 0*a1*c1))/(2. _d 0*a1)
423           IF (Tsf .GT. 0. _d 0) THEN            dTsrf(i,j) = (flxNet + k12(i,j)*(tIc1(i,j)-Tsf(i,j)))
424         &                  /(k12(i,j)-dFlxdT(i,j))
425              Tsf(i,j) = Tsf(i,j) + dTsrf(i,j)
426              IF ( Tsf(i,j) .GT. 0. _d 0 ) THEN
427  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
428              IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)             IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)
429       &        'ThSI_SOLVE4T: k,ts,t1,dTs=', k,Tsf,Tice(1),dTsf       &     'ThSI_SOLVE4T: k,ts,t1,dTs=',k,Tsf(i,j),tIc1(i,j),dTsrf(i,j)
430  #endif  #endif
431              a1 = a10 + k12             a1 = a10(i,j) + k12(i,j)
432              b1 = b10          ! note b1 = b10 - k12*Tf0  C      note: b1 = b10 - k12*Tf0
433              Tice(1) = (-b1 - SQRT(b1*b1-4. _d 0*a1*c1))/(2. _d 0*a1)             b1 = b10(i,j)
434              Tsf = 0. _d 0             tIc1(i,j) = (-b1 - SQRT(b1*b1-4. _d 0*a1*c1))/(2. _d 0*a1)
435               Tsf(i,j) = 0. _d 0
436             IF ( useBlkFlx ) THEN             IF ( useBlkFlx ) THEN
437              IF (hs.GT.3. _d -1) THEN              flxTexSW(i,j) = flx0exSW(i,j)
438                   iceornot=2              evapT(i,j) = evap0(i,j)
439              ELSE              dTsrf(i,j) = 0. _d 0
                  iceornot=1  
             ENDIF  
             IF ( useBulkForce ) THEN  
              CALL THSICE_GET_BULKF(  
      I                            iceornot, Tsf,  
      O                            flxExceptSw, df0dT, evap, dEvdT,  
      I                            i,j,bi,bj,myThid )  
             ELSEIF ( useEXF ) THEN  
              CALL THSICE_GET_EXF  (  
      I                            iceornot, Tsf,  
      O                            flxExceptSw, df0dT, evap, dEvdT,  
      I                            i,j,bi,bj,myThid )  
             ENDIF  
             dTsf = 0. _d 0  
440             ELSE             ELSE
441              flxExceptSw = flxExSW(i,j,0)              flxTexSW(i,j) = flxExSW(i,j,0)
442              dTsf = 1000.              dTsrf(i,j) = 1000.
443              df0dT = 0.              dFlxdT(i,j) = 0.
444             ENDIF             ENDIF
445             flx0 = fswdn + flxExceptSw            ENDIF
          ENDIF  
446    
447  C Check for convergence.  If no convergence, then repeat.  C--   Check for convergence.  If no convergence, then repeat.
448  C            iceFlag(i,j) = ABS(dTsrf(i,j)).GE.Terrmax
449  C Convergence test: Make sure Tsfc has converged, within prescribed error.            iterate4Tsf = iterate4Tsf .OR. iceFlag(i,j)
450  C (Energy conservation is guaranteed within machine roundoff, even  
451  C if Tsfc has not converged.)  C     Convergence test: Make sure Tsfc has converged, within prescribed error.
452  C If no convergence, then repeat.  C     (Energy conservation is guaranteed within machine roundoff, even
453    C      if Tsfc has not converged.)
454    C     If no convergence, then repeat.
455    
456  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
457           IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)            IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)
458       &     'ThSI_SOLVE4T: k,ts,t1,dTs=', k,Tsf,Tice(1),dTsf       &    'ThSI_SOLVE4T: k,ts,t1,dTs=', k,Tsf(i,j),tIc1(i,j),dTsrf(i,j)
459              IF ( useBlkFlx .AND. k.EQ.nitMaxTsf .AND. iceFlag(i,j) ) THEN
460               WRITE (6,'(A,4I4,I12,F15.9)')
461         &       ' BB: not converge: i,j,it,hi=',i,j,bi,bj,myIter,hIce(i,j)
462               WRITE (6,*) 'BB: not converge: Tsf, dTsf=',
463         &          Tsf(i,j), dTsrf(i,j)
464               WRITE (6,*) 'BB: not converge: flxNet,dFlxT=',
465         &          flxNet, dFlxdT(i,j)
466               IF (Tsf(i,j).LT.-70. _d 0) STOP
467              ENDIF
468  #endif  #endif
          IF ( useBlkFlx .AND. k.EQ.nitMaxTsf  
      &                  .AND. ABS(dTsf).GE.Terrmax ) THEN  
             WRITE (6,'(A,4I4,I12,F15.9)')  
      &                 ' BB: not converge: i,j,it,hi=',i,j,bi,bj,  
      &                                                    myIter,hi  
             WRITE (6,*) 'BB: not converge: Tsf, dTsf=', Tsf,dTsf  
             WRITE (6,*) 'BB: not converge: flx0,dfdT=',flx0,df0dT  
             IF (Tsf.LT.-70. _d 0) STOP  
          ENDIF  
469    
470             ENDIF
471            ENDDO
472           ENDDO
473    #ifndef ALLOW_AUTODIFF_TAMC
474         ENDIF         ENDIF
475    #endif
476        ENDDO        ENDDO
477  C ------ end iteration ------------  C ------ end iteration ------------
478    
479  C Compute new bottom layer temperature.  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
480    
481  #ifdef ALLOW_AUTODIFF_TAMC        DO j = jMin, jMax
482  CADJ STORE  Tice(:)      = comlev1_thsice_1, key=ikey_1         DO i = iMin, iMax
483  CADJ STORE  df0dt        = comlev1_thsice_1, key=ikey_1          IF ( iceMask(i,j).GT.0. _d 0 ) THEN
484  #endif  
485        Tice(2) = (2. _d 0*dt*k32*(Tice(1)+2. _d 0*Tf)  C--   Compute new bottom layer temperature.
486       &        + rhoi*cpIce *hi*Tice(2))            k32       = 2. _d 0*kIce  / hIce(i,j)
487       &         /(6. _d 0*dt*k32 + rhoi*cpIce *hi)            tIc2(i,j) = ( 2. _d 0*dt*k32*(tIc1(i,j)+2. _d 0*tFrz(i,j))
488         &                 + rhoi*cpIce*hIce(i,j)*tIc2(i,j))
489         &               /(6. _d 0*dt*k32 + rhoi*cpIce*hIce(i,j))
490  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
491        IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)            IF ( dBug(i,j,bi,bj) ) WRITE(6,1010)
492       &  'ThSI_SOLVE4T: k, Ts, Tice=',k,Tsf,Tice       &  'ThSI_SOLVE4T: k, Ts, Tice=',k,Tsf(i,j),tIc1(i,j),tIc2(i,j)
493              netSW   = flxAtm(i,j)
494  #endif  #endif
495    
496  C Compute final flux values at surfaces.  C--   Compute final flux values at surfaces.
497              tSrf(i,j)  = Tsf(i,j)
498        fct    = k12*(Tsf-Tice(1))            fct    = k12(i,j)*(Tsf(i,j)-tIc1(i,j))
499        flxCnB(i,j) = 4. _d 0*kIce *(Tice(2)-Tf)/hi            flxCnB(i,j) = 4. _d 0*kIce *(tIc2(i,j)-tFrz(i,j))/hIce(i,j)
500        flx0   = flx0 + df0dT*dTsf            flxNet = sHeat(i,j) + flxTexSW(i,j)
501        IF ( useBlkFlx ) THEN            flxNet = flxNet + dFlxdT(i,j)*dTsrf(i,j)
502  C--   needs to update also Evap (Tsf changes) since Latent heat has been updated            IF ( useBlkFlx ) THEN
503          evpAtm(i,j) = evap + dEvdT*dTsf  C-    needs to update also Evap (Tsf changes) since Latent heat has been updated
504        ELSE              evpAtm(i,j) = evapT(i,j) + dEvdT(i,j)*dTsrf(i,j)
505              ELSE
506  C- WARNING: Evap & +Evap*Lfresh are missing ! (but only affects Diagnostics)  C- WARNING: Evap & +Evap*Lfresh are missing ! (but only affects Diagnostics)
507          evpAtm(i,j) = 0.              evpAtm(i,j) = 0.
508        ENDIF            ENDIF
509  C-    energy flux to Atmos: use net short-wave flux at surf. and  C-    Update energy flux to Atmos with other than SW contributions;
510  C     use latent heat = Lvap (energy=0 for liq. water at 0.oC)  C     use latent heat = Lvap (energy=0 for liq. water at 0.oC)
511        flxAtm(i,j) = netSW + flxExceptSw            flxAtm(i,j) = flxAtm(i,j) + flxTexSW(i,j)
512       &                    + df0dT*dTsf + evpAtm(i,j)*Lfresh       &                + dFlxdT(i,j)*dTsrf(i,j) + evpAtm(i,j)*Lfresh
513  C-    excess of energy @ surface (used for surface melting):  C-    excess of energy @ surface (used for surface melting):
514        sHeat(i,j) = flx0 - fct            sHeat(i,j) = flxNet - fct
   
 C-    SW flux at sea-ice base left to the ocean  
       flxSW(i,j) = fswocn  
515    
516  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
517        IF ( dBug(i,j,bi,bj) ) WRITE(6,1020)            IF ( dBug(i,j,bi,bj) ) WRITE(6,1020)
518       &  'ThSI_SOLVE4T: flx0,fct,Dif,flxCnB=',       &     'ThSI_SOLVE4T: flxNet,fct,Dif,flxCnB=',
519       &                 flx0,fct,flx0-fct,flxCnB(i,j)       &                    flxNet,fct,flxNet-fct,flxCnB(i,j)
520  #endif  #endif
521    
522  C Compute new enthalpy for each layer.  C--   Compute new enthalpy for each layer.
523              qIc1(i,j) = -cpWater*Tmlt1 + cpIce *(Tmlt1-tIc1(i,j))
524        qicen(1) = -cpWater*Tmlt1 + cpIce *(Tmlt1-Tice(1))       &                + Lfresh*(1. _d 0-Tmlt1/tIc1(i,j))
525       &            + Lfresh*(1. _d 0-Tmlt1/Tice(1))            qIc2(i,j) = -cpIce *tIc2(i,j) + Lfresh
       qicen(2) = -cpIce *Tice(2) + Lfresh  
   
 C Make sure internal ice temperatures do not exceed Tmlt.  
 C (This should not happen for reasonable values of i0swFrac)  
526    
       IF (Tice(1) .GE. Tmlt1) THEN  
         WRITE (6,'(A,2I4,2I3,1P2E14.6)')  
      &   ' BBerr - Bug: IceT(1) > Tmlt',i,j,bi,bj,Tice(1),Tmlt1  
       ENDIF  
       IF (Tice(2) .GE. 0. _d 0) THEN  
        WRITE (6,'(A,2I4,2I3,1P2E14.6)')  
      &   ' BBerr - Bug: IceT(2) > 0',i,j,bi,bj,Tice(2)  
       ENDIF  
   
 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  
 C--    Update Sea-Ice state :  
           tSrf(i,j)  = Tsf  
           tIc1(i,j)  = Tice(1)  
           tic2(i,j)  = Tice(2)  
           qIc1(i,j)  = qicen(1)  
           qIc2(i,j)  = qicen(2)  
 c         dTsrf(i,j) = dTsf  
           IF ( .NOT.useBlkFlx ) dTsrf(i,j) = dTsf  
 c         sHeat(i,j) = sHeating  
 c         flxCnB(i,j)= flxCnB  
 c         flxAtm(i,j)= flxAtm  
 c         evpAtm(i,j)= evpAtm  
527  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
528    C--   Make sure internal ice temperatures do not exceed Tmlt.
529    C     (This should not happen for reasonable values of i0swFrac)
530              IF (tIc1(i,j) .GE. Tmlt1) THEN
531               WRITE (6,'(A,2I4,2I3,1P2E14.6)')
532         &     ' BBerr - Bug: IceT(1) > Tmlt',i,j,bi,bj,tIc1(i,j),Tmlt1
533              ENDIF
534              IF (tIc2(i,j) .GE. 0. _d 0) THEN
535               WRITE (6,'(A,2I4,2I3,1P2E14.6)')
536         &     ' BBerr - Bug: IceT(2) > 0',i,j,bi,bj,tIc2(i,j)
537              ENDIF
538    
539            IF ( dBug(i,j,bi,bj) ) THEN            IF ( dBug(i,j,bi,bj) ) THEN
540             WRITE(6,1020) 'ThSI_SOLV_4T: Tsf, Tice(1,2), dTsurf=',             WRITE(6,1020) 'ThSI_SOLV_4T: Tsf, Tice(1,2), dTsurf=',
541       &                              Tsf, Tice, dTsf       &           Tsf(i,j), tIc1(i,j), tIc2(i,j), dTsrf(i,j)
542             WRITE(6,1020) 'ThSI_SOLV_4T: sHeat, flxCndBt, Qice =',             WRITE(6,1020) 'ThSI_SOLV_4T: sHeat, flxCndBt, Qice =',
543       &           sHeat(i,j), flxCnB(i,j), qicen       &           sHeat(i,j), flxCnB(i,j), qIc1(i,j), qIc2(i,j)
544             WRITE(6,1020) 'ThSI_SOLV_4T: flxA, evpA, fxSW_bf,af=',             WRITE(6,1020) 'ThSI_SOLV_4T: flxA, evpA, fxSW_bf,af=',
545       &           flxAtm(i,j), evpAtm(i,j), netSW, flxSW(i,j)       &           flxAtm(i,j), evpAtm(i,j), netSW, flxSW(i,j)
546            ENDIF            ENDIF
547  #endif  #endif
548    
549          ELSE          ELSE
550            IF ( .NOT.useBlkFlx ) dTsrf(i,j) = 0. _d 0  C--     ice-free grid point:
551    c         tIc1  (i,j) = 0. _d 0
552    c         tIc2  (i,j) = 0. _d 0
553              dTsrf (i,j) = 0. _d 0
554    c         sHeat (i,j) = 0. _d 0
555    c         flxCnB(i,j) = 0. _d 0
556    c         flxAtm(i,j) = 0. _d 0
557    c         evpAtm(i,j) = 0. _d 0
558    
559          ENDIF          ENDIF
560         ENDDO         ENDDO
561        ENDDO        ENDDO
Legend:
Removed from v.1.17  
changed lines
  Added in v.1.22
  ViewVC Help
Powered by ViewVC 1.1.22