/[MITgcm]/MITgcm/pkg/thsice/thsice_calc_thickn.F
ViewVC logotype

Diff of /MITgcm/pkg/thsice/thsice_calc_thickn.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.4 by jmc, Fri Dec 17 03:44:52 2004 UTC revision 1.5 by jmc, Fri Feb 10 00:24:12 2006 UTC
# Line 41  C     compact  :: fraction of grid area Line 41  C     compact  :: fraction of grid area
41  C     hi       :: ice height  C     hi       :: ice height
42  C     hs       :: snow height  C     hs       :: snow height
43  C     Tsf      :: surface (ice or snow) temperature  C     Tsf      :: surface (ice or snow) temperature
44  C     qicen    :: ice enthalpy (J m-3)  C     qicen    :: ice enthalpy (J/kg)
45  C     qleft    :: net heat flux to ocean    [W/m2]          (> 0 downward)  C     qleft    :: net heat flux to ocean    [W/m2]          (> 0 downward)
46  C     fresh    :: Total fresh water flux to ocean [kg/m2/s] (> 0 downward)  C     fresh    :: Total fresh water flux to ocean [kg/m2/s] (> 0 downward)
47  C     fsalt    :: salt flux to ocean        [g/m2/s]        (> 0 downward)  C     fsalt    :: salt flux to ocean        [g/m2/s]        (> 0 downward)
# Line 103  C     evap           ::  evaporation ove Line 103  C     evap           ::  evaporation ove
103        _RL  dhs           ! change in snow thickness        _RL  dhs           ! change in snow thickness
104        _RL  rq            ! rho * q for a layer        _RL  rq            ! rho * q for a layer
105        _RL  rqh           ! rho * q * h for a layer        _RL  rqh           ! rho * q * h for a layer
106        _RL  qbot          ! q for new ice at bottom surf (J m-3)        _RL  qbot          ! enthalpy for new ice at bottom surf (J/kg)
107        _RL  dt            ! timestep        _RL  dt            ! timestep
108        _RL  esurp         ! surplus energy from melting (J m-2)        _RL  esurp         ! surplus energy from melting (J m-2)
109        _RL  mwater0       ! fresh water mass gained/lost (kg/m^2)        _RL  mwater0       ! fresh water mass gained/lost (kg/m^2)
# Line 134  C....................................... Line 134  C.......................................
134  C.. Compute growth and/or melting at the top and bottom surfaces.......  C.. Compute growth and/or melting at the top and bottom surfaces.......
135  C......................................................................  C......................................................................
136    
137        if (frzmlt.ge. 0. _d 0) then        IF (frzmlt.GE. 0. _d 0) THEN
138  C     !-----------------------------------------------------------------  C     !-----------------------------------------------------------------
139  C     ! freezing conditions  C     ! freezing conditions
140  C     !-----------------------------------------------------------------  C     !-----------------------------------------------------------------
141  C if higher than hihig, use all frzmlt energy to grow extra ice  C if higher than hihig, use all frzmlt energy to grow extra ice
142          if (hi.gt.hihig.and. compact.le.iceMaskmax) then          IF (hi.GT.hihig .AND. compact.LE.iceMaskmax) THEN
143            Fbot=0. _d 0            Fbot=0. _d 0
144          else          ELSE
145            Fbot=frzmlt            Fbot=frzmlt
146          endif          ENDIF
147        else        ELSE
148  C     !-----------------------------------------------------------------  C     !-----------------------------------------------------------------
149  C     ! melting conditions  C     ! melting conditions
150  C     !-----------------------------------------------------------------  C     !-----------------------------------------------------------------
151           ustar = 5. _d -2        !for no currents           ustar = 5. _d -2        !for no currents
152  C frictional velocity between ice and water  C frictional velocity between ice and water
153           ustar = sqrt(0.00536 _d 0*oceV2s)           ustar = SQRT(0.00536 _d 0*oceV2s)
154           ustar=max(5. _d -3,ustar)           ustar=max(5. _d -3,ustar)
155           cpchr =cpwater*rhosw*transcoef           cpchr =cpwater*rhosw*transcoef
156           Fbot = cpchr*(Tf-oceTs)*ustar  ! < 0           Fbot = cpchr*(Tf-oceTs)*ustar  ! < 0
157           Fbot = max(Fbot,frzmlt)    ! frzmlt < Fbot < 0           Fbot = max(Fbot,frzmlt)    ! frzmlt < Fbot < 0
158           Fbot = min(Fbot,0. _d 0)           Fbot = min(Fbot,0. _d 0)
159        endif        ENDIF
160    
161  C  mass of fresh water and salt initially present in ice  C  mass of fresh water and salt initially present in ice
162        mwater0 = rhos*hs + rhoi*hi        mwater0 = rhos*hs + rhoi*hi
# Line 169  C---+----1----+----2----+----3----+----4 Line 169  C---+----1----+----2----+----3----+----4
169    
170  C Compute energy available for melting/growth.  C Compute energy available for melting/growth.
171    
172        if (hi.lt.himin0) then        IF (hi.LT.himin0) THEN
173  C below a certain height, all energy goes to changing ice extent  C below a certain height, all energy goes to changing ice extent
174         frace=1. _d 0         frace=1. _d 0
175        else        ELSE
176         frace=frac_energy         frace=frac_energy
177        endif        ENDIF
178        if (hi.gt.hihig) then        IF (hi.GT.hihig) THEN
179  C above certain height only melt from top  C above certain height only melt from top
180         frace=0. _d 0         frace=0. _d 0
181        else        ELSE
182         frace=frac_energy         frace=frac_energy
183        endif        ENDIF
184  C force this when no ice fractionation  C force this when no ice fractionation
185        if (frac_energy.eq.0. _d 0) frace=0. _d 0        IF (frac_energy.EQ.0. _d 0) frace=0. _d 0
186    
187  c     if (Tsf .eq. 0. _d 0 .and. sHeating.gt.0. _d 0) then  c     IF (Tsf .EQ. 0. _d 0 .AND. sHeating.GT.0. _d 0) THEN
188        if ( sHeating.gt.0. _d 0 ) then        IF ( sHeating.GT.0. _d 0 ) THEN
189            etop = (1. _d 0-frace)*sHeating * dt            etop = (1. _d 0-frace)*sHeating * dt
190            etope = frace*sHeating * dt            etope = frace*sHeating * dt
191        else        ELSE
192            etop =  0. _d 0            etop =  0. _d 0
193            etope = 0. _d 0            etope = 0. _d 0
194  C jmc: found few cases where Tsf=0 & sHeating < 0 : add this line to conserv energy:  C jmc: found few cases where Tsf=0 & sHeating < 0 : add this line to conserv energy:
195            esurp = sHeating * dt            esurp = sHeating * dt
196        endif        ENDIF
197  C--   flux at the base of sea-ice:  C--   flux at the base of sea-ice:
198  C     conduction H.flx= flxCnB (+ =down); oceanic turbulent H.flx= Fbot (+ =down).  C     conduction H.flx= flxCnB (+ =down); oceanic turbulent H.flx= Fbot (+ =down).
199  C-    ==> energy available(+ => melt)= (flxCnB-Fbot)*dt  C-    ==> energy available(+ => melt)= (flxCnB-Fbot)*dt
200  c     if (frzmlt.lt.0. _d 0) then  c     IF (frzmlt.LT.0. _d 0) THEN
201  c         ebot = (1. _d 0-frace)*(flxCnB-Fbot) * dt  c         ebot = (1. _d 0-frace)*(flxCnB-Fbot) * dt
202  c         ebote = frace*(flxCnB-Fbot) * dt  c         ebote = frace*(flxCnB-Fbot) * dt
203  c     else  c     ELSE
204  c         ebot = (flxCnB-Fbot) * dt  c         ebot = (flxCnB-Fbot) * dt
205  c         ebote = 0. _d 0  c         ebote = 0. _d 0
206  c     endif  c     ENDIF
207  C- original formulation(above): Loose energy when flxCnB < Fbot < 0  C- original formulation(above): Loose energy when flxCnB < Fbot < 0
208        ebot = (flxCnB-Fbot) * dt        ebot = (flxCnB-Fbot) * dt
209        if (ebot.gt.0. _d 0) then        IF (ebot.GT.0. _d 0) THEN
210           ebote = frace*ebot           ebote = frace*ebot
211           ebot  = ebot-ebote           ebot  = ebot-ebote
212        else        ELSE
213           ebote = 0. _d 0           ebote = 0. _d 0
214        endif        ENDIF
215        IF (dBug) WRITE(6,1020) 'ThSI_CALC_TH: etop,etope,ebot,ebote=',        IF (dBug) WRITE(6,1020) 'ThSI_CALC_TH: etop,etope,ebot,ebote=',
216       &            etop,etope,ebot,ebote       &            etop,etope,ebot,ebote
217    
# Line 221  C If they do, then eliminate the layer. Line 221  C If they do, then eliminate the layer.
221  C for reasonable values of i0.)  C for reasonable values of i0.)
222    
223        hlyr = hi / rnlyr        hlyr = hi / rnlyr
224        do k = 1, nlyr        DO k = 1, nlyr
225           hnew(k) = hlyr           hnew(k) = hlyr
226        enddo        ENDDO
227    
228  C Top melt: snow, then ice.  C Top melt: snow, then ice.
229    
230        if (etop .gt. 0. _d 0) then        IF (etop .GT. 0. _d 0) THEN
231           if (hs. gt. 0. _d 0) then           IF (hs. gt. 0. _d 0) THEN
232              rq =  rhos * qsnow              rq =  rhos * qsnow
233              rqh = rq * hs              rqh = rq * hs
234              if (etop .lt. rqh) then              IF (etop .LT. rqh) THEN
235                 hs = hs - etop/rq                 hs = hs - etop/rq
236                 etop = 0. _d 0                 etop = 0. _d 0
237              else              ELSE
238                 etop = etop - rqh                 etop = etop - rqh
239                 hs = 0. _d 0                 hs = 0. _d 0
240              endif              ENDIF
241           endif           ENDIF
242                                        
243           do k = 1, nlyr           DO k = 1, nlyr
244              if (etop .gt. 0. _d 0) then              IF (etop .GT. 0. _d 0) THEN
245                 rq =  rhoi * qicen(k)                 rq =  rhoi * qicen(k)
246                 rqh = rq * hnew(k)                 rqh = rq * hnew(k)
247                 if (etop .lt. rqh) then                 IF (etop .LT. rqh) THEN
248                    hnew(k) = hnew(k) - etop / rq                    hnew(k) = hnew(k) - etop / rq
249                    etop = 0. _d 0                    etop = 0. _d 0
250                 else                 ELSE
251                    etop = etop - rqh                    etop = etop - rqh
252                    hnew(k) = 0. _d 0                    hnew(k) = 0. _d 0
253                 endif                 ENDIF
254              endif              ENDIF
255           enddo           ENDDO
256        else        ELSE
257          etop=0. _d 0          etop=0. _d 0
258        endif        ENDIF
259  C If ice is gone and melting energy remains  C If ice is gone and melting energy remains
260  c     if (etop .gt. 0. _d 0) then  c     IF (etop .GT. 0. _d 0) THEN
261  c        write (6,*)  'QQ All ice melts from top  ', i,j  c        WRITE (6,*)  'QQ All ice melts from top  ', i,j
262  c        hi=0. _d 0  c        hi=0. _d 0
263  c        go to 200  c        go to 200
264  c     endif  c     ENDIF
265    
266    
267  C Bottom melt/growth.  C Bottom melt/growth.
268    
269        if (ebot .lt. 0. _d 0) then        IF (ebot .LT. 0. _d 0) THEN
270  C Compute enthalpy of new ice growing at bottom surface.  C Compute enthalpy of new ice growing at bottom surface.
271           qbot =  -cpice *Tf + Lfresh           qbot =  -cpice *Tf + Lfresh
272           dhi = -ebot / (qbot * rhoi)           dhi = -ebot / (qbot * rhoi)
# Line 274  C Compute enthalpy of new ice growing at Line 274  C Compute enthalpy of new ice growing at
274           k = nlyr           k = nlyr
275           qicen(k) = (hnew(k)*qicen(k)+dhi*qbot) / (hnew(k)+dhi)           qicen(k) = (hnew(k)*qicen(k)+dhi*qbot) / (hnew(k)+dhi)
276           hnew(k) = hnew(k) + dhi           hnew(k) = hnew(k) + dhi
277        else                      ELSE              
278           do k = nlyr, 1, -1           DO k = nlyr, 1, -1
279              if (ebot.gt.0. _d 0 .and. hnew(k).gt.0. _d 0) then              IF (ebot.GT.0. _d 0 .AND. hnew(k).GT.0. _d 0) THEN
280                 rq =  rhoi * qicen(k)                 rq =  rhoi * qicen(k)
281                 rqh = rq * hnew(k)                 rqh = rq * hnew(k)
282                 if (ebot .lt. rqh) then                 IF (ebot .LT. rqh) THEN
283                    hnew(k) = hnew(k) - ebot / rq                    hnew(k) = hnew(k) - ebot / rq
284                    ebot = 0. _d 0                    ebot = 0. _d 0
285                 else                 ELSE
286                    ebot = ebot - rqh                    ebot = ebot - rqh
287                    hnew(k) = 0. _d 0                    hnew(k) = 0. _d 0
288                 endif                 ENDIF
289              endif              ENDIF
290           enddo           ENDDO
291    
292  C If ice melts completely and snow is left, remove the snow with  C If ice melts completely and snow is left, remove the snow with
293  C energy from the mixed layer  C energy from the mixed layer
294    
295           if (ebot.gt.0. _d 0 .and. hs.gt.0. _d 0) then           IF (ebot.GT.0. _d 0 .AND. hs.GT.0. _d 0) THEN
296              rq =  rhos * qsnow              rq =  rhos * qsnow
297              rqh = rq * hs              rqh = rq * hs
298              if (ebot .lt. rqh) then              IF (ebot .LT. rqh) THEN
299                 hs = hs - ebot / rq                 hs = hs - ebot / rq
300                 ebot = 0. _d 0                 ebot = 0. _d 0
301              else              ELSE
302                 ebot = ebot - rqh                 ebot = ebot - rqh
303                 hs = 0. _d 0                 hs = 0. _d 0
304              endif              ENDIF
305           endif           ENDIF
306  c        if (ebot .gt. 0. _d 0) then  c        IF (ebot .GT. 0. _d 0) THEN
307  c           IF (dBug) WRITE(6,*) 'All ice (& snow) melts from bottom'  c           IF (dBug) WRITE(6,*) 'All ice (& snow) melts from bottom'
308  c           hi=0. _d 0  c           hi=0. _d 0
309  c           go to 200  c           go to 200
310  c        endif  c        ENDIF
311        endif        ENDIF
312    
313  C Compute new total ice thickness.  C Compute new total ice thickness.
314        hi = 0. _d 0        hi = 0. _d 0
315        do k = 1, nlyr        DO k = 1, nlyr
316           hi = hi + hnew(k)           hi = hi + hnew(k)
317        enddo        ENDDO
318        IF (dBug) WRITE(6,1020) 'ThSI_CALC_TH:   etop, ebot, hi, hs =',        IF (dBug) WRITE(6,1020) 'ThSI_CALC_TH:   etop, ebot, hi, hs =',
319       &                                         etop, ebot, hi, hs           &                                         etop, ebot, hi, hs    
320    
321  C If hi < himin, melt the ice.  C If hi < himin, melt the ice.
322        if ( hi.lt.himin .AND. (hi+hs).gt.0. _d 0 ) then        IF ( hi.LT.himin .AND. (hi+hs).GT.0. _d 0 ) THEN
323           esurp = esurp - rhos*qsnow*hs           esurp = esurp - rhos*qsnow*hs
324           do k = 1, nlyr           DO k = 1, nlyr
325              esurp = esurp - rhoi*qicen(k)*hnew(k)              esurp = esurp - rhoi*qicen(k)*hnew(k)
326           enddo           ENDDO
327           hi = 0. _d 0           hi = 0. _d 0
328           hs = 0. _d 0           hs = 0. _d 0
329           Tsf=0. _d 0           Tsf=0. _d 0
330           compact = 0. _d 0           compact = 0. _d 0
331           do k = 1, nlyr           DO k = 1, nlyr
332             qicen(k) = 0. _d 0             qicen(k) = 0. _d 0
333           enddo           ENDDO
334           IF (dBug) WRITE(6,1020) 'ThSI_CALC_TH: -1 : esurp=',esurp           IF (dBug) WRITE(6,1020) 'ThSI_CALC_TH: -1 : esurp=',esurp
335        endif        ENDIF
336    
337  C--   do a mass-budget of sea-ice to compute "fresh" = the fresh-water flux  C--   do a mass-budget of sea-ice to compute "fresh" = the fresh-water flux
338  C     that is returned to the ocean ; needs to be done before snow/evap  C     that is returned to the ocean ; needs to be done before snow/evap
# Line 340  C     that is returned to the ocean ; ne Line 340  C     that is returned to the ocean ; ne
340    
341  C- note : was not supposed to modify snowPr in THSICE_CALC_TH ;  C- note : was not supposed to modify snowPr in THSICE_CALC_TH ;
342  C         but to reproduce old results, reset to zero if Tsf >= 0  C         but to reproduce old results, reset to zero if Tsf >= 0
343  c     IF ( Tsf.ge.0. _d 0 ) snowPr = 0.  c     IF ( Tsf.GE.0. _d 0 ) snowPr = 0.
344    
345        IF ( hi .LE. 0. _d 0 ) THEN        IF ( hi .LE. 0. _d 0 ) THEN
346  C-    return snow to the ocean (account for Latent heat of freezing)  C-    return  snow to the ocean (account for Latent heat of freezing)
347          fresh = fresh + snowPr          fresh = fresh + snowPr
348          qleft = qleft - snowPr*Lfresh          qleft = qleft - snowPr*Lfresh
349          GOTO 200  
350        ENDIF        ELSE
351    C-    else: hi > 0
352    
353  C Let it snow  C Let it snow
354    
# Line 355  C Let it snow Line 356  C Let it snow
356    
357  C If ice evap is used to sublimate surface snow/ice or  C If ice evap is used to sublimate surface snow/ice or
358  C if no ice pass on to ocean  C if no ice pass on to ocean
359        if (hs.gt.0. _d 0) then         IF (hs.GT.0. _d 0) THEN
360          if (evap/rhos *dt.gt.hs) then           IF (evap/rhos *dt.GT.hs) THEN
361             evap=evap-hs*rhos/dt             evap=evap-hs*rhos/dt
362             hs=0. _d 0             hs=0. _d 0
363          else           ELSE
364             hs = hs - evap/rhos *dt             hs = hs - evap/rhos *dt
365             evap=0. _d 0             evap=0. _d 0
366          endif           ENDIF
367        endif         ENDIF
368        if (hi.gt.0. _d 0.and.evap.gt.0. _d 0) then         IF (hi.GT.0. _d 0.AND.evap.GT.0. _d 0) THEN
369           do k = 1, nlyr           DO k = 1, nlyr
370              if (evap .gt. 0. _d 0) then              IF (evap .GT. 0. _d 0) THEN
371  C-- original scheme, does not care about ice temp.  C-- original scheme, does not care about ice temp.
372  C-  this can produce small error (< 1.W/m2) in the Energy budget  C-  this can produce small error (< 1.W/m2) in the Energy budget
373  c              if (evap/rhoi *dt.gt.hnew(k)) then  c              IF (evap/rhoi *dt.GT.hnew(k)) THEN
374  c                evap=evap-hnew(k)*rhoi/dt  c                evap=evap-hnew(k)*rhoi/dt
375  c                hnew(k)=0. _d 0  c                hnew(k)=0. _d 0
376  c              else  c              ELSE
377  c                hnew(k) = hnew(k) - evap/rhoi *dt  c                hnew(k) = hnew(k) - evap/rhoi *dt
378  c                evap=0. _d 0  c                evap=0. _d 0
379  c              endif  c              ENDIF
380  C-- modified scheme. taking into account Ice enthalpy  C-- modified scheme. taking into account Ice enthalpy
381                 dhi = evap/rhoi*dt                 dhi = evap/rhoi*dt
382                 if (dhi.ge.hnew(k)) then                 IF (dhi.GE.hnew(k)) THEN
383                   evap=evap-hnew(k)*rhoi/dt                   evap=evap-hnew(k)*rhoi/dt
384                   esurp = esurp - hnew(k)*rhoi*(qicen(k)-Lfresh)                   esurp = esurp - hnew(k)*rhoi*(qicen(k)-Lfresh)
385                   hnew(k)=0. _d 0                   hnew(k)=0. _d 0
386                 else                 ELSE
387                   hq = hnew(k)*qicen(k)-dhi*Lfresh                   hq = hnew(k)*qicen(k)-dhi*Lfresh
388                   hnew(k) = hnew(k) - dhi                   hnew(k) = hnew(k) - dhi
389                   qicen(k)=hq/hnew(k)                   qicen(k)=hq/hnew(k)
390                   evap=0. _d 0                   evap=0. _d 0
391                 endif                 ENDIF
392  C-------  C-------
393              endif              ENDIF
394           enddo           ENDDO
395        endif         ENDIF
396  c     if (evap .gt. 0. _d 0) then  c     IF (evap .GT. 0. _d 0) THEN
397  c           write (6,*)  'BB All ice sublimates', i,j  c           WRITE (6,*)  'BB All ice sublimates', i,j
398  c           hi=0. _d 0  c           hi=0. _d 0
399  c           go to 200  c           go to 200
400  c     endif  c     ENDIF
401    
402  C Compute new total ice thickness.  C Compute new total ice thickness.
403    
404        hi = 0. _d 0         hi = 0. _d 0
405        do k = 1, nlyr         DO k = 1, nlyr
406           hi = hi + hnew(k)           hi = hi + hnew(k)
407        enddo         ENDDO
408    
409  C If hi < himin, melt the ice.  C If hi < himin, melt the ice.
410        if ( hi.gt.0. _d 0 .AND. hi.lt.himin ) then         IF ( hi.GT.0. _d 0 .AND. hi.LT.himin ) THEN
411           fresh = fresh + (rhos*hs + rhoi*hi)/dt           fresh = fresh + (rhos*hs + rhoi*hi)/dt
412           esurp = esurp - rhos*qsnow*hs           esurp = esurp - rhos*qsnow*hs
413           do k = 1, nlyr           DO k = 1, nlyr
414              esurp = esurp - rhoi*qicen(k)*hnew(k)              esurp = esurp - rhoi*qicen(k)*hnew(k)
415           enddo           ENDDO
416           hi = 0. _d 0           hi = 0. _d 0
417           hs = 0. _d 0           hs = 0. _d 0
418           Tsf=0. _d 0           Tsf=0. _d 0
419           compact = 0. _d 0           compact = 0. _d 0
420           do k = 1, nlyr           DO k = 1, nlyr
421             qicen(k) = 0. _d 0             qicen(k) = 0. _d 0
422           enddo           ENDDO
423           IF (dBug) WRITE(6,1020) 'ThSI_CALC_TH: -2 : esurp,fresh=',           IF (dBug) WRITE(6,1020) 'ThSI_CALC_TH: -2 : esurp,fresh=',
424       &                   esurp, fresh       &                   esurp, fresh
425        endif         ENDIF
426        IF ( hi .le. 0. _d 0 ) GOTO 200  
427    C-    else hi > 0: end
428          ENDIF
429    
430          IF ( hi .GT. 0. _d 0 ) THEN
431    
432  C If there is enough snow to lower the ice/snow interface below  C If there is enough snow to lower the ice/snow interface below
433  C freeboard, convert enough snow to ice to bring the interface back  C freeboard, convert enough snow to ice to bring the interface back
434  C to sea-level.  Adjust enthalpy of top ice layer accordingly.  C to sea-level.  Adjust enthalpy of top ice layer accordingly.
435    
436        if ( hs .gt. hi*rhoiw/rhos ) then         IF ( hs .GT. hi*rhoiw/rhos ) THEN
437  cBB               write (6,*)  'Freeboard adjusts'  cBB               WRITE (6,*)  'Freeboard adjusts'
438           dhi = (hs * rhos - hi * rhoiw) / rhosw           dhi = (hs * rhos - hi * rhoiw) / rhosw
439           dhs = dhi * rhoi / rhos           dhs = dhi * rhoi / rhos
440           rqh = rhoi*qicen(1)*hnew(1) + rhos*qsnow*dhs           rqh = rhoi*qicen(1)*hnew(1) + rhos*qsnow*dhs
# Line 437  cBB               write (6,*)  'Freeboar Line 442  cBB               write (6,*)  'Freeboar
442           qicen(1) = rqh / (rhoi*hnew(1))           qicen(1) = rqh / (rhoi*hnew(1))
443           hi = hi + dhi           hi = hi + dhi
444           hs = hs - dhs           hs = hs - dhs
445        end if         ENDIF
446    
447    
448  C limit ice height  C limit ice height
449  C- NOTE: this part does not conserve Energy ;  C- NOTE: this part does not conserve Energy ;
450  C        but surplus of fresh water and salt are taken into account.  C        but surplus of fresh water and salt are taken into account.
451        if (hi.gt.hiMax) then         IF (hi.GT.hiMax) THEN
452  cBB      print*,'BBerr, hi>hiMax',i,j,hi  cBB      print*,'BBerr, hi>hiMax',i,j,hi
453           chi=hi-hiMax           chi=hi-hiMax
454           do k=1,nlyr           DO k=1,nlyr
455              hnew(k)=hnew(k)-chi/2. _d 0              hnew(k)=hnew(k)-chi/2. _d 0
456           enddo           ENDDO
457           fresh = fresh + chi*rhoi/dt           fresh = fresh + chi*rhoi/dt
458        endif         ENDIF
459        if (hs.gt.hsMax) then         IF (hs.GT.hsMax) THEN
460  c        print*,'BBerr, hs>hsMax',i,j,hs  c        print*,'BBerr, hs>hsMax',i,j,hs
461           chs=hs-hsMax           chs=hs-hsMax
462           hs=hsMax           hs=hsMax
463           fresh = fresh + chs*rhos/dt           fresh = fresh + chs*rhos/dt
464        endif         ENDIF
465    
466  C Compute new total ice thickness.  C Compute new total ice thickness.
467    
468        hi = 0. _d 0         hi = 0. _d 0
469        do k = 1, nlyr         DO k = 1, nlyr
470           hi = hi + hnew(k)           hi = hi + hnew(k)
471        enddo         ENDDO
472    
473        IF (dBug) WRITE(6,1020) 'ThSI_CALC_TH: b-Winton: hnew, qice =',         IF (dBug) WRITE(6,1020) 'ThSI_CALC_TH: b-Winton: hnew, qice =',
474       &                                                 hnew, qicen       &                                                  hnew, qicen
475    
476         hlyr = hi/rnlyr         hlyr = hi/rnlyr
477         CALL THSICE_RESHAPE_LAYERS(         CALL THSICE_RESHAPE_LAYERS(
478       U                            qicen,       U                            qicen,
479       I                            hlyr, hnew, myThid )       I                            hlyr, hnew, myThid )
480    
481        IF (dBug) WRITE(6,1020) 'ThSI_CALC_TH: compact,hi, qtot, hs =',         IF (dBug) WRITE(6,1020) 'ThSI_CALC_TH: compact,hi, qtot, hs =',
482       &                  compact,hi,(qicen(1)+qicen(2))*0.5, hs       &                  compact,hi,(qicen(1)+qicen(2))*0.5, hs
483    
484  200   continue  C-    if hi > 0 : end
485          ENDIF
486     200  CONTINUE
487    
488  C-  Compute surplus energy left over from melting.  C-  Compute surplus energy left over from melting.
489    
490        if (hi.le.0. _d 0) compact=0. _d 0        IF (hi.LE.0. _d 0) compact=0. _d 0
491    
492  C.. heat fluxes left over for ocean  C.. heat fluxes left over for ocean
493         qleft = qleft + (Fbot+(esurp+etop+ebot)/dt)         qleft = qleft + (Fbot+(esurp+etop+ebot)/dt)
# Line 515  C      and it can only be reduced by lat Line 522  C      and it can only be reduced by lat
522    
523  C calculate extent changes  C calculate extent changes
524        extend=etope+ebote        extend=etope+ebote
525        if (compact.gt.0. _d 0.and.extend.gt.0. _d 0) then        IF (compact.GT.0. _d 0.AND.extend.GT.0. _d 0) THEN
526           rq =  rhoi * 0.5 _d 0*(qicen(1)+qicen(2))           rq =  rhoi * 0.5 _d 0*(qicen(1)+qicen(2))
527           rs =  rhos * qsnow           rs =  rhos * qsnow
528           rqh = rq * hi + rs * hs           rqh = rq * hi + rs * hs
529           freshe=(rhos*hs+rhoi*hi)/dt           freshe=(rhos*hs+rhoi*hi)/dt
530           salte=(rhoi*hi*saltice)/dt           salte=(rhoi*hi*saltice)/dt
531           if (extend .lt. rqh) then           IF (extend .LT. rqh) THEN
532             compact=(1. _d 0-extend/rqh)*compact             compact=(1. _d 0-extend/rqh)*compact
533             fresh=fresh+extend/rqh*freshe             fresh=fresh+extend/rqh*freshe
534             fsalt=fsalt+extend/rqh*salte             fsalt=fsalt+extend/rqh*salte
535           else           ELSE
536             compact=0. _d 0             compact=0. _d 0
537             hi=0. _d 0             hi=0. _d 0
538             hs=0. _d 0             hs=0. _d 0
539             qleft=qleft+(extend-rqh)/dt             qleft=qleft+(extend-rqh)/dt
540             fresh=fresh+freshe             fresh=fresh+freshe
541             fsalt=fsalt+salte             fsalt=fsalt+salte
542           endif           ENDIF
543        elseif (extend.gt.0. _d 0) then        ELSEIF (extend.GT.0. _d 0) THEN
544           qleft=qleft+extend/dt           qleft=qleft+extend/dt
545        endif        ENDIF
546    
547  #endif  /* ALLOW_THSICE */  #endif  /* ALLOW_THSICE */
548    
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.22