File: C:\NOAA\NEMS_11731\src\atmos\phys\sfc_sice.f

1     !-----------------------------------
2           subroutine sfc_sice                                               &
3     !...................................
4     !  ---  inputs:
5          &     ( im, km, ps, u1, v1, t1, q1, delt,                          &
6          &       sfcemis, dlwflx, sfcnsw, sfcdsw, srflag,                   &
7          &       cm, ch, prsl1, prslki, slimsk, ddvel,                      &
8          &       flag_iter, mom4ice, lsm, lprnt,ipr,                        &
9     !  ---  input/outputs:
10          &       hice, fice, tice, weasd, tskin, tprcp, stc, ep,            &
11     !  ---  outputs:
12          &       snwdph, qsurf, snowmt, gflux, cmm, chh, evap, hflx         &
13          &     )
14     
15     ! ===================================================================== !
16     !  description:                                                         !
17     !                                                                       !
18     !  usage:                                                               !
19     !                                                                       !
20     !    call sfc_sice                                                      !
21     !       inputs:                                                         !
22     !          ( im, km, ps, u1, v1, t1, q1, delt,                          !
23     !            sfcemis, dlwflx, sfcnsw, sfcdsw, srflag,                   !
24     !            cm, ch, prsl1, prslki, slimsk, ddvel,                      !
25     !            flag_iter, mom4ice, lsm,                                   !
26     !       input/outputs:                                                  !
27     !            hice, fice, tice, weasd, tskin, tprcp, stc, ep,            !
28     !       outputs:                                                        !
29     !            snwdph, qsurf, snowmt, gflux, cmm, chh, evap, hflx )       !
30     !                                                                       !
31     !  subprogram called:  ice3lay.                                         !
32     !                                                                       !
33     !  program history log:                                                 !
34     !         2005  --  xingren wu created  from original progtm and added  !
35     !                     two-layer ice model                               !
36     !         200x  -- sarah lu    added flag_iter                          !
37     !    oct  2006  -- h. wei      added cmm and chh to output              !
38     !         2007  -- x. wu modified for mom4 coupling (i.e. mom4ice)      !
39     !         2007  -- s. moorthi micellaneous changes                      !
40     !    may  2009  -- y.-t. hou   modified to include surface emissivity   !
41     !                     effect on lw radiation. replaced the confusing    !
42     !                     slrad with sfc net sw sfcnsw (dn-up). reformatted !
43     !                     the code and add program documentation block.     !
44     !    sep  2009 -- s. moorthi removed rcl, changed pressure units and    !
45     !                     further optimized                                 !
46     !                                                                       !
47     !                                                                       !
48     !  ====================  defination of variables  ====================  !
49     !                                                                       !
50     !  inputs:                                                       size   !
51     !     im, km   - integer, horiz dimension and num of soil layers   1    !
52     !     ps       - real, surface pressure                            im   !
53     !     u1, v1   - real, u/v component of surface layer wind         im   !
54     !     t1       - real, surface layer mean temperature ( k )        im   !
55     !     q1       - real, surface layer mean specific humidity        im   !
56     !     delt     - real, time interval (second)                      1    !
57     !     sfcemis  - real, sfc lw emissivity ( fraction )              im   !
58     !     dlwflx   - real, total sky sfc downward lw flux ( w/m**2 )   im   !
59     !     sfcnsw   - real, total sky sfc netsw flx into ground(w/m**2) im   !
60     !     sfcdsw   - real, total sky sfc downward sw flux ( w/m**2 )   im   !
61     !     srflag   - real, snow/rain flag for precipitation            im   !
62     !     cm       - real, surface exchange coeff for momentum (m/s)   im   !
63     !     ch       - real, surface exchange coeff heat & moisture(m/s) im   !
64     !     prsl1    - real, surface layer mean pressure                 im   !
65     !     prslki   - real,                                             im   !
66     !     slimsk   - real, sea/land/ice mask (=0/1/2)                  im   !
67     !     ddvel    - real,                                             im   !
68     !     flag_iter- logical,                                          im   !
69     !     mom4ice  - logical,                                          im   !
70     !     lsm      - integer, flag for land surface model scheme       1    !
71     !                =0: use osu scheme; =1: use noah scheme                !
72     !                                                                       !
73     !  input/outputs:                                                       !
74     !     hice     - real, sea-ice thickness                           im   !
75     !     fice     - real, sea-ice concentration                       im   !
76     !     tice     - real, sea-ice surface temperature                 im   !
77     !     weasd    - real, water equivalent accumulated snow depth (mm)im   !
78     !     tskin    - real, ground surface skin temperature ( k )       im   !
79     !     tprcp    - real, total precipitation                         im   !
80     !     stc      - real, soil temp (k)                              im,km !
81     !     ep       - real, potential evaporation                       im   !
82     !                                                                       !
83     !  outputs:                                                             !
84     !     snwdph   - real, water equivalent snow depth (mm)            im   !
85     !     qsurf    - real, specific humidity at sfc                    im   !
86     !     snowmt   - real, snow melt (m)                               im   !
87     !     gflux    - real, soil heat flux (w/m**2)                     im   !
88     !     cmm      - real,                                             im   !
89     !     chh      - real,                                             im   !
90     !     evap     - real, evaperation from latent heat flux           im   !
91     !     hflx     - real, sensible heat flux                          im   !
92     !                                                                       !
93     ! ===================================================================== !
94     !
95           use machine , only : kind_phys
96           use funcphys, only : fpvs
97           use physcons, only : sbc => con_sbc, hvap => con_hvap,            &
98          &                     tgice => con_tice, cp => con_cp,             &
99          &                     eps => con_eps, epsm1 => con_epsm1,          &
100          &                     grav => con_g, rvrdm1 => con_fvirt,          &
101          &                     t0c => con_t0c, rd => con_rd
102     !
103           implicit none
104     !
105     !  ---  constant parameters:
106           integer,              parameter :: kmi   = 2        ! 2-layer of ice
107           real(kind=kind_phys), parameter :: cpinv = 1.0/cp
108           real(kind=kind_phys), parameter :: hvapi = 1.0/hvap
109           real(kind=kind_phys), parameter :: elocp = hvap/cp
110           real(kind=kind_phys), parameter :: himax = 8.0      ! maximum ice thickness allowed
111           real(kind=kind_phys), parameter :: himin = 0.1      ! minimum ice thickness required
112           real(kind=kind_phys), parameter :: hsmax = 2.0      ! maximum snow depth allowed
113           real(kind=kind_phys), parameter :: timin = 173.0    ! minimum temperature allowed for snow/ice
114           real(kind=kind_phys), parameter :: albfw = 0.06     ! albedo for lead
115           real(kind=kind_phys), parameter :: dsi   = 1.0/0.33
116     
117     !  ---  inputs:
118           integer, intent(in) :: im, km, lsm, ipr
119           logical, intent(in) :: lprnt
120     
121           real (kind=kind_phys), dimension(im), intent(in) :: ps, u1, v1,   &
122          &       t1, q1, sfcemis, dlwflx, sfcnsw, sfcdsw, srflag, cm, ch,   &
123          &       prsl1, prslki, slimsk, ddvel
124     
125           real (kind=kind_phys), intent(in) :: delt
126     
127           logical, intent(in) :: flag_iter(im), mom4ice
128     
129     !  ---  input/outputs:
130           real (kind=kind_phys), dimension(im), intent(inout) :: hice,      &
131          &       fice, tice, weasd, tskin, tprcp, ep
132     
133           real (kind=kind_phys), dimension(im,km), intent(inout) :: stc
134     
135     !  ---  outputs:
136           real (kind=kind_phys), dimension(im), intent(out) :: snwdph,      &
137          &       qsurf, snowmt, gflux, cmm, chh, evap, hflx
138     
139     !  ---  locals:
140           real (kind=kind_phys), dimension(im) :: ffw, evapi, evapw,        &
141          &       hflxi, hflxw, sneti, snetw, qssi, qssw, hfd, hfi, hfw,     &
142          &       focn, snof, hi_save, hs_save, psurf, q0, qs1, rch, rho,    &
143          &       snowd, theta1, tv1, ps1, wind
144     
145           real (kind=kind_phys) :: cimin, t12, t14, tem, stsice(im,kmi)
146     
147           integer :: i, k
148      
149           logical :: flag(im)
150     !
151     !===> ...  begin here
152     !
153     !  --- ...  set minimum ice concentration required
154     
155           if (mom4ice) then
156              cimin = 0.15          ! mom4ice and mask
157           else
158              cimin = 0.50          ! gfs only
159           endif
160     
161     !  --- ...  set flag for sea-ice
162     
163           do i = 1, im
164              flag(i) = (slimsk(i)>=2.0) .and. flag_iter(i)
165           enddo
166     
167           if (mom4ice) then
168             do i = 1, im
169               if (flag(i)) then
170                 hi_save(i) = hice(i)
171                 hs_save(i) = weasd(i) * 0.001
172               endif
173             enddo
174           endif
175     !
176           do i = 1, im
177             if (flag_iter(i) .and. slimsk(i)<1.5) then
178               hice(i) = 0.0
179               fice(i) = 0.0
180             endif
181           enddo
182     
183     !  --- ...  snow-rain detection
184     
185           if (.not. mom4ice .and. lsm > 0) then
186             do i = 1, im
187               if (flag(i)) then
188                 if (srflag(i) == 1.0) then
189                   ep(i) = 0.0
190                   weasd(i) = weasd(i) + 1.e3*tprcp(i)
191                   tprcp(i)  = 0.0
192                 endif
193               endif
194             enddo
195           endif
196     
197     !  --- ...  initialize variables. all units are supposedly m.k.s. unless specifie
198     !           psurf is in pascals, wind is wind speed, theta1 is adiabatic surface
199     !           temp from level 1, rho is density, qs1 is sat. hum. at level1 and qss
200     !           is sat. hum. at surface
201     !           convert slrad to the civilized unit from langley minute-1 k-4
202     
203           do i = 1, im
204             if (flag(i)) then
205               psurf(i) = 1000.0 * ps(i)
206               ps1(i)   = 1000.0 * prsl1(i)
207     
208     !         dlwflx has been given a negative sign for downward longwave
209     !         sfcnsw is the net shortwave flux (direction: dn-up)
210     
211               wind(i)   = sqrt(u1(i)*u1(i) + v1(i)*v1(i))                   &
212          &              + max(0.0, min(ddvel(i), 30.0))
213               wind(i)   = max(wind(i), 1.0)
214     
215               q0(i)     = max(q1(i), 1.0e-8)
216     !         tsurf(i)  = tskin(i)
217               theta1(i) = t1(i) * prslki(i)
218               tv1(i)    = t1(i) * (1.0 + rvrdm1*q0(i))
219               rho(i)    = prsl1(i) / (rd*tv1(i))
220               qs1(i)    = fpvs(t1(i))
221               qs1(i)    = eps*qs1(i) / (prsl1(i) + epsm1*qs1(i))
222               qs1(i)    = max(qs1(i), 1.e-8)
223               q0 (i)    = min(qs1(i), q0(i))
224     
225               ffw(i)    = 1.0 - fice(i)
226               if (fice(i) < cimin) then
227                 print *,'warning: ice fraction is low:', fice(i)
228                 fice(i) = cimin
229                 ffw (i) = 1.0 - fice(i)
230                 tice(i) = tgice
231                 tskin(i)= tgice
232                 print *,'fix ice fraction: reset it to:', fice(i)
233               endif
234     
235               qssi(i) = fpvs(tice(i))
236               qssi(i) = eps*qssi(i) / (ps(i) + epsm1*qssi(i))
237               qssw(i) = fpvs(tgice)
238               qssw(i) = eps*qssw(i) / (ps(i) + epsm1*qssw(i))
239     
240     !  --- ...  snow depth in water equivalent is converted from mm to m unit
241     
242               if (mom4ice) then
243                 snowd(i) = weasd(i) * 0.001 / fice(i)
244               else
245                 snowd(i) = weasd(i) * 0.001
246               endif
247     !         flagsnw(i) = .false.
248     
249     !  --- ...  when snow depth is less than 1 mm, a patchy snow is assumed and
250     !           soil is allowed to interact with the atmosphere.
251     !           we should eventually move to a linear combination of soil and
252     !           snow under the condition of patchy snow.
253             endif
254           enddo
255     
256     
257           do i = 1, im
258             if (flag(i)) then
259     
260     !  --- ...  rcp = rho cp ch v
261     
262               cmm(i) = cm(i)  * wind(i)
263               chh(i) = rho(i) * ch(i) * wind(i)
264               rch(i) = chh(i) * cp
265     
266     !  --- ...  sensible and latent heat flux over open water & sea ice
267     
268               evapi(i) = elocp * rch(i) * (qssi(i) - q0(i))
269               evapw(i) = elocp * rch(i) * (qssw(i) - q0(i))
270     !         evap(i)  = fice(i)*evapi(i) + ffw(i)*evapw(i)
271             endif
272           enddo
273     
274     !  --- ...  update sea ice temperature
275     
276           do k = 1, kmi
277             do i = 1, im
278               if (flag(i)) then
279                 stsice(i,k) = stc(i,k)
280               endif
281             enddo
282           enddo
283     
284           if (lprnt) write(0,*)' tice=',tice(ipr)
285     
286           do i = 1, im
287             if (flag(i)) then
288               snetw(i) = sfcdsw(i) * (1.0 - albfw)
289               snetw(i) = min(3.0*sfcnsw(i)/(1.0+2.0*ffw(i)), snetw(i))
290               sneti(i) = (sfcnsw(i) - ffw(i)*snetw(i)) / fice(i)
291     
292               t12 = tice(i) * tice(i)
293               t14 = t12 * t12
294     
295     !  --- ...  hfi = net non-solar and upir heat flux @ ice surface
296     
297               hfi(i) = -dlwflx(i) + sfcemis(i)*sbc*t14 + evapi(i)           &
298          &           + rch(i)*(tice(i) - theta1(i))
299               hfd(i) = 4.0*sfcemis(i)*sbc*tice(i)*t12                       &
300          &           + (1.0 + elocp*eps*hvap*qs1(i)/(rd*t12)) * rch(i)
301     
302               t12 = tgice * tgice
303               t14 = t12 * t12
304     
305     !  --- ...  hfw = net heat flux @ water surface (within ice)
306     
307     !         hfw(i) = -dlwflx(i) + sfcemis(i)*sbc*t14 + evapw(i)           &
308     !    &           + rch(i)*(tgice - theta1(i)) - snetw(i)
309     
310               focn(i) = 2.0     ! heat flux from ocean - should be from ocn model
311               snof(i) = 0.0     ! snowfall rate - snow accumulates in gbphys
312     
313               hice(i) = max( min( hice(i), himax ), himin )
314               snowd(i) = min( snowd(i), hsmax )
315     
316               if (snowd(i) > (2.0*hice(i))) then
317                 print *, 'warning: too much snow :',snowd(i)
318                 snowd(i) = 2.0 * hice(i)
319                 print *,'fix: decrease snow depth to:',snowd(i)
320               endif
321             endif
322           enddo
323     
324           if (lprnt) write(0,*)' tice2=',tice(ipr)
325           call ice3lay
326     !  ---  inputs:                                                         !
327     !    &     ( im, kmi, fice, flag, hfi, hfd, sneti, focn, delt,          !
328     !  ---  outputs:                                                        !
329     !    &       snowd, hice, stsice, tice, snof, snowmt, gflux )           !
330     
331           if (lprnt) write(0,*)' tice3=',tice(ipr)
332           if (mom4ice) then
333             do i = 1, im
334               if (flag(i)) then
335                 hice(i)  = hi_save(i)
336                 snowd(i) = hs_save(i)
337               endif
338             enddo
339           endif
340     
341           do i = 1, im
342             if (flag(i)) then
343               if (tice(i) < timin) then
344                 print *,'warning: snow/ice temperature is too low:',tice(i)
345          &,' i=',i
346                 tice(i) = timin
347                 print *,'fix snow/ice temperature: reset it to:',tice(i)
348               endif
349     
350               if (stsice(i,1) < timin) then
351                 print *,'warning: layer 1 ice temp is too low:',stsice(i,1)
352          &,' i=',i
353                 stsice(i,1) = timin
354                 print *,'fix layer 1 ice temp: reset it to:',stsice(i,1)
355               endif
356     
357               if (stsice(i,2) < timin) then
358                 print *,'warning: layer 2 ice temp is too low:',stsice(i,2)
359                 stsice(i,2) = timin
360                 print *,'fix layer 2 ice temp: reset it to:',stsice(i,2)
361               endif
362     
363               tskin(i) = tice(i)*fice(i) + tgice*ffw(i)
364             endif
365           enddo
366     
367           do k = 1, kmi
368             do i = 1, im
369               if (flag(i)) then
370                 stc(i,k) = min(stsice(i,k), t0c)
371               endif
372             enddo
373           enddo
374     
375     !  --- ...  calculate sensible heat flux (& evap over sea ice)
376     
377           do i = 1, im
378             if (flag(i)) then
379               hflxi(i) = rch(i) * (tice(i) - theta1(i))
380               hflxw(i) = rch(i) * (tgice - theta1(i))
381               hflx(i)  = fice(i)*hflxi(i) + ffw(i)*hflxw(i)
382               evap(i)  = fice(i)*evapi(i) + ffw(i)*evapw(i)
383             endif
384           enddo
385     !
386     !  --- ...  the rest of the output
387     
388           do i = 1, im
389             if (flag(i)) then
390               qsurf(i) = q1(i) + evap(i) / (elocp*rch(i))
391     
392     !  --- ...  convert snow depth back to mm of water equivalent
393     
394               weasd(i)  = snowd(i) * 1000.0
395               snwdph(i) = weasd(i) * dsi             ! snow depth in mm
396             endif
397           enddo
398     
399           do i = 1, im
400             if (flag(i)) then
401               tem     = 1.0 / rho(i)
402               hflx(i) = hflx(i) * tem * cpinv
403               evap(i) = evap(i) * tem * hvapi
404             endif
405           enddo
406     !
407           return
408     
409     ! =================
410           contains
411     ! =================
412     
413     
414     !-----------------------------------
415           subroutine ice3lay
416     !...................................
417     !  ---  inputs:
418     !    &     ( im, kmi, fice, flag, hfi, hfd, sneti, focn, delt,          &
419     !  ---  input/outputs:
420     !    &       snowd, hice, stsice, tice, snof,                           &
421     !  ---  outputs:
422     !    &       snowmt, gflux                                              &
423     !    &     )
424     
425     !**************************************************************************
426     !                                                                         *
427     !            three-layer sea ice vertical thermodynamics                  *
428     !                                                                         *
429     ! based on:  m. winton, "a reformulated three-layer sea ice model",       *
430     ! journal of atmospheric and oceanic technology, 2000                     *
431     !                                                                         *
432     !                                                                         *
433     !        -> +---------+ <- tice - diagnostic surface temperature ( <= 0c )*
434     !       /   |         |                                                   *
435     !   snowd   |  snow   | <- 0-heat capacity snow layer                     *
436     !       \   |         |                                                   *
437     !        => +---------+                                                   *
438     !       /   |         |                                                   *
439     !      /    |         | <- t1 - upper 1/2 ice temperature; this layer has *
440     !     /     |         |         a variable (t/s dependent) heat capacity  *
441     !   hice    |...ice...|                                                   *
442     !     \     |         |                                                   *
443     !      \    |         | <- t2 - lower 1/2 ice temp. (fixed heat capacity) *
444     !       \   |         |                                                   *
445     !        -> +---------+ <- base of ice fixed at seawater freezing temp.   *
446     !                                                                         *
447     !  =====================  defination of variables  =====================  !
448     !                                                                         !
449     !  inputs:                                                         size   !
450     !     im, kmi  - integer, horiz dimension and num of ice layers      1    !
451     !     fice     - real, sea-ice concentration                         im   !
452     !     flag     - logical, ice mask flag                              1    !
453     !     hfi      - real, net non-solar and heat flux @ surface(w/m^2)  im   !
454     !     hfd      - real, heat flux derivatice @ sfc (w/m^2/deg-c)      im   !
455     !     sneti    - real, net solar incoming at top  (w/m^2)            im   !
456     !     focn     - real, heat flux from ocean    (w/m^2)               im   !
457     !     delt     - real, timestep                (sec)                 1    !
458     !                                                                         !
459     !  input/outputs:                                                         !
460     !     snowd    - real, surface pressure                              im   !
461     !     hice     - real, sea-ice thickness                             im   !
462     !     stsice   - real, temp @ midpt of ice levels  (deg c)          im,kmi!     
463     !     tice     - real, surface temperature     (deg c)               im   !
464     !     snof     - real, snowfall rate           (m/sec)               im   !
465     !                                                                         !
466     !  outputs:                                                               !
467     !     snowmt   - real, snow melt during delt   (m)                   im   !
468     !     gflux    - real, conductive heat flux    (w/m^2)               im   !
469     !                                                                         !
470     !  locals:                                                                !
471     !     hdi      - real, ice-water interface     (m)                   im   !
472     !     hsni     - real, snow-ice                (m)                   im   !
473     !                                                                         !
474     ! ======================================================================= !
475     !
476     
477     !  ---  constant parameters: (properties of ice, snow, and seawater)
478           real (kind=kind_phys), parameter :: ds   = 330.0    ! snow (ov sea ice) density (kg/m^3)
479           real (kind=kind_phys), parameter :: dw   =1000.0    ! fresh water density  (kg/m^3)
480           real (kind=kind_phys), parameter :: t0c  =273.15    ! freezing temp of fresh ice (k)
481           real (kind=kind_phys), parameter :: ks   = 0.31     ! conductivity of snow   (w/mk)
482           real (kind=kind_phys), parameter :: i0   = 0.3      ! ice surface penetrating solar fraction
483           real (kind=kind_phys), parameter :: ki   = 2.03     ! conductivity of ice  (w/mk)
484           real (kind=kind_phys), parameter :: di   = 917.0    ! density of ice   (kg/m^3)
485           real (kind=kind_phys), parameter :: ci   = 2054.0   ! heat capacity of fresh ice (j/kg/k)
486           real (kind=kind_phys), parameter :: li   = 3.34e5   ! latent heat of fusion (j/kg-ice)
487           real (kind=kind_phys), parameter :: si   = 1.0      ! salinity of sea ice
488           real (kind=kind_phys), parameter :: mu   = 0.054    ! relates freezing temp to salinity
489           real (kind=kind_phys), parameter :: tfi  = -mu*si   ! sea ice freezing temp = -mu*salinity
490           real (kind=kind_phys), parameter :: tfw  = -1.8     ! tfw - seawater freezing temp (c)
491           real (kind=kind_phys), parameter :: tfi0 = tfi-0.0001
492           real (kind=kind_phys), parameter :: dici = di*ci 
493           real (kind=kind_phys), parameter :: dili = di*li 
494           real (kind=kind_phys), parameter :: dsli = ds*li 
495           real (kind=kind_phys), parameter :: ki4  = ki*4.0
496     
497     !  ---  inputs:
498     !     integer, intent(in) :: im, kmi
499     
500     !     real (kind=kind_phys), dimension(im), intent(in) :: fice, hfi,    &
501     !    &       hfd, sneti, focn
502     
503     !     real (kind=kind_phys), intent(in) :: delt
504     
505     !     logical, dimension(im), intent(in) :: flag
506     
507     !  ---  input/outputs:
508     !     real (kind=kind_phys), dimension(im), intent(inout) :: snowd,     &
509     !    &       hice, tice, snof
510     
511     !     real (kind=kind_phys), dimension(im,kmi), intent(inout) :: stsice
512     
513     !  ---  outputs:
514     !     real (kind=kind_phys), dimension(im), intent(out) :: snowmt,      &
515     !    &       gflux
516     
517     !  ---  locals:
518           real (kind=kind_phys), dimension(im) :: hdi, hsni, a, b, ip,      &
519          &      a1, b1, c1, a10, b10, k12, k32, h1, h2, dh, f1, tsf,        &
520          &      tmelt, bmelt
521     
522           real (kind=kind_phys) :: dt2, dt4, dt6
523     
524           integer :: i
525     !
526     !===> ...  begin here
527     !
528           dt2 = 2.0 * delt
529           dt4 = 4.0 * delt
530           dt6 = 6.0 * delt
531     
532           do i = 1, im
533             if (flag(i)) then
534               snowd(i) = snowd(i) * dw / ds
535               hdi(i) = (ds*snowd(i) + di*hice(i)) / dw
536     
537               if (hice(i) < hdi(i)) then
538                 snowd(i) = snowd(i) + hice(i) - hdi(i)
539                 hsni (i) = (hdi(i) - hice(i)) * ds / di
540                 hice (i) = hice(i) + hsni(i)
541               endif
542     
543               snof(i) = snof(i) * dw / ds
544               tice(i) = tice(i) - t0c
545               stsice(i,1) = min(stsice(i,1)-t0c, tfi0)     ! degc
546               stsice(i,2) = min(stsice(i,2)-t0c, tfi0)     ! degc
547     
548               ip(i) = i0 * sneti(i)         ! ip +v (in winton ip=-i0*sneti as sol -v)
549               if (snowd(i) > 0.0) then
550                 tsf(i) = 0.0
551                 ip (i) = 0.0
552               else
553                 tsf(i) = tfi
554                 ip (i) = i0 * sneti(i)      ! ip +v here (in winton ip=-i0*sneti)
555               endif
556               tice(i) = min(tice(i), tsf(i))
557     
558     !  --- ...  compute ice temperature
559     
560               b(i) = hfd(i)
561               a(i) = hfi(i) - sneti(i) + ip(i) - tice(i)*b(i)  ! +v sol input here
562               k12(i) = ki4*ks / (ks*hice(i) + ki4*snowd(i))
563               k32(i) = 2.0 * ki / hice(i)
564     
565               a10(i) = dici*hice(i)/dt2 + k32(i)*(dt4*k32(i) + dici*hice(i))&
566          &           / (dt6*k32(i) + dici*hice(i))
567               b10(i) = -di*hice(i) * (ci*stsice(i,1) + li*tfi/stsice(i,1))  &
568          &           / dt2 - ip(i)                                          &
569          &           - k32(i)*(dt4*k32(i)*tfw + dici*hice(i)*stsice(i,2))   &
570          &           / (dt6*k32(i) + dici*hice(i))
571     
572               a1(i) = a10(i) + k12(i)*b(i) / (k12(i) + b(i))
573               b1(i) = b10(i) + a(i)*k12(i) / (k12(i) + b(i))
574               c1(i) = dili*tfi / dt2*hice(i)
575     
576               stsice(i,1) = -(sqrt(b1(i)*b1(i) - 4.0*a1(i)*c1(i))           &
577          &                + b1(i))/(2.0*a1(i))
578               tice(i) = (k12(i)*stsice(i,1) - a(i)) / (k12(i) + b(i))
579     
580               if (tice(i) > tsf(i)) then
581                 a1(i) = a10(i) + k12(i)
582                 b1(i) = b10(i) - k12(i)*tsf(i)
583                 stsice(i,1) = -(sqrt(b1(i)*b1(i) - 4.0*a1(i)*c1(i))         &
584          &                  + b1(i)) / (2.0*a1(i))
585                 tice(i) = tsf(i)
586                 tmelt(i) = (k12(i)*(stsice(i,1) - tsf(i))                   &
587          &               - (a(i) + b(i)*tsf(i))) * delt
588               else
589                 tmelt(i) = 0.0
590                 snowd(i) = snowd(i) + snof(i)*delt
591               endif
592     
593               stsice(i,2) = (dt2*k32(i)*(stsice(i,1) + 2.0*tfw)             &
594          &                + dici*hice(i)*stsice(i,2))                       &
595          &                / (dt6*k32(i) + dici*hice(i))
596     
597               bmelt(i) = (focn(i) + ki4*(stsice(i,2) - tfw)/hice(i)) * delt
598     
599     !  --- ...  resize the ice ...
600     
601               h1(i) = 0.5 * hice(i)
602               h2(i) = 0.5 * hice(i)
603     
604     !  --- ...  top ...
605     
606               if (tmelt(i) <= snowd(i)*dsli) then
607                 snowmt(i) = tmelt(i) / dsli
608                 snowd (i) = snowd(i) - tmelt(i)/dsli
609               else
610                 snowmt(i) = snowd(i)
611                 h1(i) = h1(i) - (tmelt(i) - snowd(i)*dsli)                  &
612          &            / (di * (ci - li/stsice(i,1)) * (tfi - stsice(i,1)))
613                 snowd(i) = 0.0
614               endif
615     
616     !  --- ...  and bottom
617     
618               if (bmelt(i) < 0.0) then
619                 dh(i) = -bmelt(i) / (dili + dici*(tfi - tfw))
620                 stsice(i,2) = (h2(i)*stsice(i,2) + dh(i)*tfw)               &
621          &                  / (h2(i) + dh(i))
622                 h2(i) = h2(i) + dh(i)
623               else
624                 h2(i) = h2(i) - bmelt(i) / (dili + dici*(tfi - stsice(i,2)))
625               endif
626     
627     !  --- ...  if ice remains, even up 2 layers, else, pass negative energy back in snow
628     
629               hice(i) = h1(i) + h2(i)
630     
631               if (hice(i) > 0.0) then
632                 if (h1(i) > 0.5*hice(i)) then
633                   f1(i) = 1.0 - 2.0*h2(i) / hice(i)
634                   stsice(i,2) = f1(i)                                       &
635          &                    * (stsice(i,1) + li*tfi/(ci*stsice(i,1)))     &
636          &                    + (1.0 - f1(i))*stsice(i,2)
637     
638                   if (stsice(i,2) > tfi) then
639                     hice(i) = hice(i) - h2(i)*ci*(stsice(i,2) - tfi)        &
640          &                  / (li*delt)
641                     stsice(i,2) = tfi
642                   endif
643                 else
644                   f1(i) = 2.0 * h1(i) / hice(i)
645                   stsice(i,1) = f1(i)                                       &
646          &                    * (stsice(i,1) + li*tfi/(ci*stsice(i,1)))     &
647          &                    + (1.0 - f1(i))*stsice(i,2)
648                   stsice(i,1) = (stsice(i,1) - sqrt(stsice(i,1)*stsice(i,1) &
649          &                    - 4.0*tfi*li/ci)) / 2.0
650                 endif
651     
652                 k12(i) = ki4*ks / (ks*hice(i) + ki4*snowd(i))
653                 gflux(i) = k12(i) * (stsice(i,1) - tice(i))
654               else
655                 snowd(i) = snowd(i) + (h1(i)*(ci*(stsice(i,1) - tfi)        &
656          &               - li*(1.0 - tfi/stsice(i,1)))                      &
657          &               + h2(i)*(ci*(stsice(i,2) - tfi) - li)) / li
658     
659                 hice(i) = max(0.0, snowd(i)*ds/di)
660                 snowd(i) = 0.0
661                 stsice(i,1) = tfw
662                 stsice(i,2) = tfw
663                 gflux(i) = 0.0
664               endif   ! end if_hice_block
665     
666               gflux(i) = fice(i) * gflux(i)
667               snowmt(i) = snowmt(i) * ds / dw
668               snowd(i) = snowd(i) * ds / dw
669               tice(i) = tice(i) + t0c
670               stsice(i,1) = stsice(i,1) + t0c
671               stsice(i,2) = stsice(i,2) + t0c
672             endif   ! end if_flag_block
673           enddo   ! end do_i_loop
674     
675           return
676     !...................................
677           end subroutine ice3lay
678     !-----------------------------------
679     
680     ! =========================== !
681     !     end contain programs    !
682     ! =========================== !
683     
684     !...................................
685           end subroutine sfc_sice
686     !-----------------------------------
687