My Project
ice_therm_vertical.f90
Go to the documentation of this file.
1 
2 
3 
4 
5 
6 
7 
8 
9 
10 
11 
12 !/===========================================================================/
13 ! CVS VERSION INFORMATION
14 ! $Id$
15 ! $Name$
16 ! $Revision$
17 !/===========================================================================/
18 
19 !=========================================================================
20 !BOP
21 !
22 ! !MODULE: ice_therm_vertical - thermo calculations before call to coupler
23 !
24 ! !DESCRIPTION:
25 !
26 ! Update ice and snow internal temperatures and compute
27 ! thermodynamic growth rates and atmospheric fluxes.
28 !
29 ! See Bitz, C.M., and W.H. Lipscomb, 1999:
30 ! An energy-conserving thermodynamic model of sea ice,
31 ! J. Geophys. Res., 104, 15,669-15,677.
32 !
33 ! NOTE: The thermodynamic calculation is split in two for load balancing.
34 ! First ice\_therm\_vertical computes vertical growth rates and coupler
35 ! fluxes. Then ice\_therm\_itd does thermodynamic calculations not
36 ! needed for coupling.
37 !
38 ! !REVISION HISTORY:
39 !
40 ! authors: William H. Lipscomb (LANL) and C.M. Bitz (UW)
41 !
42 ! Vectorized by Clifford Chen (Fujitsu) and William H. Lipscomb (LANL)
43 !
44 ! !INTERFACE:
45 !
46  module ice_therm_vertical
47 !
48 ! !USES:
49 !
50  use ice_model_size
51  use ice_kinds_mod
52  use ice_domain
53  use ice_fileunits
54  use ice_constants
55  use ice_calendar
56  use ice_grid
57  use ice_state
58  use ice_flux
59  use ice_itd
60  use mod_utils
61 
62 
63 
64 ! use ice_diagnostics, only: print_state
65 ! use ice_diagnostics ! print_state
66 
67 !
68 !EOP
69 !
70  implicit none
71  save
72 
73  real (kind=dbl_kind), parameter :: &
74  ferrmax = 1.0e-3_dbl_kind & ! max allowed energy flux error (W m-2)
75  ! rec&ommend ferrmax < 0.01 W m-2
76  , hsnomin = 1.0e-6_dbl_kind ! min thickness for which Tsno computed (m)
77 
78  real (kind=dbl_kind) :: &
79  salin(nilyr+1) &! salinity (ppt)
80  , tmlt(nilyr+1) &! melting temp, -mu * salinity
81  , ustar_scale ! scaling for i ce-ocean heat flux
82 
83 ! real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,ncat) ::&
84  real (kind=dbl_kind), dimension (:,:,:),allocatable,save ::&
85  hicen_old ! old ice thick ness (m), used by ice_therm_itd
86 
87 ! real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi) ::&
88  real (kind=dbl_kind), dimension (:,:) ,allocatable,save ::&
89  rside ! fraction of i ce that melts laterally
90 
91  logical (kind=log_kind) :: &
92  l_brine ! if true, treat brine pocket effects
93 
94  character (char_len), private :: stoplabel
95 
96 !=======================================================================
97 
98  contains
99 
100 !=======================================================================
101 !BOP
102 !
103 ! !ROUTINE: init_thermo_vertical - initialize salinity and melting temp
104 !
105 ! !DESCRIPTION:
106 !
107 ! Initialize the vertical profile of ice salinity and melting temperature.
108 !
109 ! !REVISION HISTORY:
110 !
111 ! authors: C. M. Bitz, UW
112 ! William H. Lipscomb, LANL
113 !
114 ! !INTERFACE:
115 !
116  subroutine init_thermo_vertical
117 !
118 ! !USES:
119 !
120  use ice_itd
121 !
122 ! !INPUT/OUTPUT PARAMETERS:
123 !
124 !EOP
125 !
126  real (kind=dbl_kind), parameter ::&
127  nsal = 0.407_dbl_kind&
128  , msal = 0.573_dbl_kind&
129  , saltmax = 3.2_dbl_kind &! max salinity at ice base (ppt)
130  , min_salin = 0.1_dbl_kind ! threshold for brine pocket treatment
131 
132  integer (kind=int_kind) :: k ! ice layer index
133  real (kind=dbl_kind) :: zn ! normalized ice thickness
134 
135  !-----------------------------------------------------------------
136  ! Determine l_brine based on saltmax.
137  ! Thermodynamic solver will not converge if l_brine is true and
138  ! saltmax is close to zero.
139  !-----------------------------------------------------------------
140 
141  if (saltmax > min_salin) then
142  l_brine = .true.
143  else
144  l_brine = .false.
145  endif
146 
147  ! salinity and melting temperature profile
148  if (l_brine) then
149  do k = 1, nilyr
150  zn = (real(k,kind=dbl_kind)-p5) / real(nilyr,kind=dbl_kind)
151  salin(k)=(saltmax/c2i)*(c1i-cos(pi*zn**(nsal/(msal+zn))))
152 ! salin(k)=saltmax ! for isosaline ice
153  enddo
154  salin(nilyr+1) = saltmax
155  do k = 1, nilyr+1
156  tmlt(k) = -salin(k)*depresst
157  enddo
158  else
159  do k = 1, nilyr+1
160  salin(k) = c0i
161  tmlt(k) = c0i
162  enddo
163  endif
164  ustar_scale = c1i ! for nonzero currents
165 
166  end subroutine init_thermo_vertical
167 
168 !=======================================================================
169 !BOP
170 !
171 ! !ROUTINE: thermo_vertical - driver for pre-coupler thermodynamics
172 !
173 ! !DESCRIPTION:
174 !
175 ! Driver for updating ice and snow internal temperatures and
176 ! computing thermodynamic growth rates and atmospheric fluxes.
177 !
178 ! NOTE: The wind stress is computed here for later use.
179 !
180 ! !REVISION HISTORY:
181 !
182 ! authors: William H. Lipscomb, LANL
183 ! C. M. Bitz, UW
184 !
185 ! !INTERFACE:
186 !
187  subroutine thermo_vertical
188 !
189 ! !USES:
190 !
191  use ice_atmo
192 ! use ice_timers
193 ! use ice_ocean
194 ! ggao
195  use ice_work, only: worka, workb
196 
197 !
198 ! !INPUT/OUTPUT PARAMETERS:
199 !
200 !EOP
201 !
202  integer (kind=int_kind) :: &
203  i, j &! horizontal indices
204  , ij &! horizontal index, combines i and j loops
205  , ni &! thickness category index
206  , k &! ice layer index
207  , icells ! number of cells with aicen > puny
208 
209  integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)) :: &
210  indxi, indxj ! compressed indices for cells with aicen > puny
211 
212  real (kind=dbl_kind) :: &
213  dhi & ! change in ice thickness
214  , dhs ! change in snow thickness
215 
216 ! 2D coupler variables (computed for each category, then aggregated)
217 
218  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi) :: &
219  strxn & ! air/ice zonal strss, (N/m^2)
220  , stryn & ! air/ice merdnl strss, (N/m^2)
221  , fsensn & ! surface downward sensible heat (W/m^2)
222  , flatn & ! surface downward latent heat (W/m^2)
223  , fswabsn & ! shortwave absorbed by ice (W/m^2)
224  , flwoutn & ! upward LW at surface (W/m^2)
225  , evapn & ! flux of vapor, atmos to ice (kg m-2 s-1)
226  , trefn & ! air tmp reference level (K)
227  , qrefn & ! air sp hum reference level (kg/kg)
228  , freshn & ! flux of water, ice to ocean (kg/m^2/s)
229  , fsaltn & ! flux of salt, ice to ocean (kg/m^2/s)
230  , fhnetn & ! fbot corrected for leftover energy (W/m^2)
231  , fswthrun ! SW through ice to ocean (W/m^2)
232 
233 ! 2D state variables (thickness, temperature, enthalpy)
234 
235  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi) :: &
236  hin & ! ice thickness (m)
237  , hsn & ! snow thickness (m)
238  , hlyr & ! ice layer thickness
239  , hin_init & ! initial value of hin
240  , hsn_init & ! initial value of hsn
241  , hsn_new & ! thickness of new snow (m)
242  , qsn & ! snow enthalpy, qsn < 0 (J m-3)
243  , tsn & ! internal snow temperature (deg C)
244  , tsf ! ice/snow top surface temp, same as Tsfcn (deg C)
245 
246  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nilyr) :: &
247  hnew & ! new ice layer thicknesses (m)
248  , qin & ! ice layer enthalpy, qin < 0 (J m-3)
249  , tin ! internal ice layer temperatures
250 
251 ! other 2D flux and energy variables
252 
253  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi) :: &
254  tbot & ! ice bottom surface temperature (deg C)
255  , fbot & ! ice-ocean heat flux at bottom surface (W/m^2)
256  , fsurf & ! net flux to top surface, not including fcondtop
257  , fcondtop & ! downward cond flux at top surface (W m-2)
258  , fcondbot & ! downward cond flux at bottom surface (W m-2)
259  , fswint & ! SW absorbed in ice interior, below surface (W m-2)
260  , einit & ! initial energy of melting (J m-2)
261  , efinal & ! final energy of melting (J m-2)
262  , mvap ! ice/snow mass sublimated/condensed (kg m-2)
263 
264 ! call ice_timer_start(4) ! column model
265 ! call ice_timer_start(5) ! thermodynamics
266 
267  !-----------------------------------------------------------------
268  ! Initialize coupler variables sent to the atmosphere.
269  !
270  ! We cannot initialize the ocean coupler fluxes (fresh, fsalt,
271  ! fhnet, and fswthru) because these may have changed during the
272  ! previous time step after the call to_coupler. The ocean
273  ! coupler fluxes are initialized by init_flux_ocn after
274  ! the call to_coupler.
275  !-----------------------------------------------------------------
276  call init_flux_atm
277 
278  !-----------------------------------------------------------------
279  ! Initialize diagnostic variables for history files.
280  !-----------------------------------------------------------------
281  call init_diagnostics
282 
283  !-----------------------------------------------------------------
284  ! Adjust frzmlt to account for ice-ocean heat fluxes since last
285  ! call to coupler.
286  ! Compute lateral and bottom heat fluxes.
287  !-----------------------------------------------------------------
288  call frzmlt_bottom_lateral (tbot, fbot, rside)
289 
290  !-----------------------------------------------------------------
291  ! Vertical thermodynamics: growth rates, updated ice state, and
292  ! coupler variables
293  !-----------------------------------------------------------------
294 
295  do ni = 1, ncat
296 
297  !-----------------------------------------------------------------
298  ! Save initial ice thickness (used later for ITD remapping)
299  !-----------------------------------------------------------------
300 
301  do j = jlo, jhi
302  do i = ilo, ihi
303  if (aicen(i,j,ni) > puny) then
304  hicen_old(i,j,ni) = vicen(i,j,ni) / aicen(i,j,ni)
305  else
306  hicen_old(i,j,ni) = c0i
307  endif
308  enddo ! i
309  enddo ! j
310 
311 
312  !-----------------------------------------------------------------
313  ! Initialize the single-category versions of coupler fluxes,
314  ! along with other thermodynamic arrays.
315  !-----------------------------------------------------------------
316 
317  call init_thermo_vars (strxn, stryn, trefn, qrefn, &
318  fsensn, flatn, fswabsn, flwoutn, &
319  evapn, freshn, fsaltn, fhnetn, &
320  fswthrun, fsurf, fcondtop, fcondbot, &
321  fswint, einit, efinal, mvap)
322 
323  !-----------------------------------------------------------------
324  ! prep for air-to-ice heat, momentum, radiative and water fluxes
325  !-----------------------------------------------------------------
326 
327  call atmo_boundary_layer (ni, 'ice', tsfcn(ilo:ihi,jlo:jhi,ni),&
328  strxn, stryn, trefn, qrefn, worka, workb)
329 
330  !-----------------------------------------------------------------
331  ! Identify cells with nonzero ice area
332  !-----------------------------------------------------------------
333 
334  icells = 0
335  do j = jlo, jhi
336  do i = ilo, ihi
337  if (aicen(i,j,ni) > puny) then
338  icells = icells + 1
339  indxi(icells) = i
340  indxj(icells) = j
341  endif
342  enddo ! i
343  enddo ! j
344 
345  !-----------------------------------------------------------------
346  ! Compute variables needed for vertical thermo calculation
347  !-----------------------------------------------------------------
348 
349  call init_vertical_profile &
350  (ni, icells, indxi, indxj, &
351  hin, hlyr, hsn, hin_init, hsn_init, &
352  qin, tin, qsn, tsn, tsf, einit)
353 
354  !-----------------------------------------------------------------
355  ! Compute new surface temperature and internal ice and snow
356  ! temperatures
357  !-----------------------------------------------------------------
358 
359  call temperature_changes &
360  (ni, icells, indxi, indxj, &
361  hlyr, hsn, qin, tin, qsn, tsn, &
362  tsf, tbot, fbot, fsensn, flatn, fswabsn, &
363  flwoutn, fswthrun, fhnetn, fsurf, fcondtop, fcondbot, &
364  fswint, einit)
365 
366  !-----------------------------------------------------------------
367  ! Compute growth and/or melting at the top and bottom surfaces.
368  ! Repartition ice into equal-thickness layers, conserving energy.
369  !-----------------------------------------------------------------
370 
371  call thickness_changes &
372  (ni, icells, indxi, indxj, &
373  hin, hlyr, hsn, qin, qsn, &
374  mvap, tbot, fbot, flatn, fhnetn, &
375  fsurf, fcondtop, fcondbot, efinal)
376 
377  !-----------------------------------------------------------------
378  ! Check for energy conservation by comparing the change in energy
379  ! to the net energy input
380  !-----------------------------------------------------------------
381 
382  call conservation_check_vthermo &
383  (ni, icells, indxi, indxj, &
384  fsurf, flatn, fhnetn, fswint, einit, efinal)
385 
386  !-----------------------------------------------------------------
387  ! Let it snow
388  !-----------------------------------------------------------------
389 
390  call add_new_snow &
391  (ni, icells, indxi, indxj, &
392  hsn, hsn_new, qsn, tsf)
393 
394  !-----------------------------------------------------------------
395  ! Compute fluxes of water and salt from ice to ocean
396  !-----------------------------------------------------------------
397 
398  do ij = 1, icells
399  i = indxi(ij)
400  j = indxj(ij)
401 
402  dhi = hin(i,j) - hin_init(i,j)
403  dhs = hsn(i,j) - hsn_init(i,j)
404 
405 ! evapn(i,j) = mvap(i,j)/dt ! mvap < 0 => sublimation,
406  evapn(i,j) = mvap(i,j)/dtice ! mvap < 0 => sublimation,
407  ! mvap > 0 => condensation
408  freshn(i,j) = evapn(i,j) - &
409  (rhoi*dhi + rhos*(dhs-hsn_new(i,j))) / dtice
410  fsaltn(i,j) = -rhoi*dhi*ice_ref_salinity*p001/dtice
411 ! (rhoi*dhi + rhos*(dhs-hsn_new(i,j))) / dt
412 ! fsaltn(i,j) = -rhoi*dhi*ice_ref_salinity*p001/dt
413 
414  enddo ! ij
415 
416  !-----------------------------------------------------------------
417  ! Increment area-weighted fluxes.
418  ! Note: Must not change aicen before calling this subroutine.
419  !-----------------------------------------------------------------
420 
421  call merge_fluxes (ni, &
422  strxn, stryn, fsensn, flatn, fswabsn, &
423  flwoutn, evapn, trefn, qrefn, &
424  freshn, fsaltn, fhnetn, fswthrun)
425 
426  !-----------------------------------------------------------------
427  ! Given the vertical thermo state variables (hin, hsn, qin,
428  ! qsn, compute the new ice state variables (vicen, vsnon,
429  ! eicen, esnon).
430  !-----------------------------------------------------------------
431 
432  call update_state_vthermo &
433  (ni, icells, indxi, indxj, &
434  hin, hsn, qin, qsn, tsf)
435 
436  enddo ! ncat
437 
438  !-----------------------------------------------------------------
439  ! Update mixed layer with heat from ice
440  !-----------------------------------------------------------------
441 
442 ! if (oceanmixed_ice) then
443 ! do j=jlo,jhi
444 ! do i=ilo,ihi
445 ! if (hmix(i,j) > c0i) sst(i,j) = sst(i,j) &
446 ! + (fhnet(i,j)+fswthru(i,j))*dt/(cp_ocn*rhow*hmix(i,j))
447 ! enddo
448 ! enddo
449 ! endif
450 
451 ! call ice_timer_stop(5) ! thermodynamics
452 ! call ice_timer_stop(4) ! column model
453 ! ggao
454 
455  end subroutine thermo_vertical
456 
457 !=======================================================================
458 !BOP
459 !
460 ! !ROUTINE: frzmlt_bottom_lateral - bottom and lateral heat fluxes
461 !
462 ! !DESCRIPTION:
463 !
464 ! Adjust frzmlt to account for changes in fhnet since from\_coupler.
465 ! Compute heat flux to bottom surface.
466 ! Compute fraction of ice that melts laterally.
467 !
468 ! !REVISION HISTORY:
469 !
470 ! authors C. M. Bitz, UW
471 ! William H. Lipscomb, LANL
472 ! Elizabeth C. Hunke, LANL
473 !
474 ! !INTERFACE:
475 !
476  subroutine frzmlt_bottom_lateral (Tbot, fbot, rside)
477 !
478 ! !USES:
479 !
480 ! !INPUT/OUTPUT PARAMETERS:
481 !
482  real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out) :: &
483  tbot & ! ice bottom surface temperature (deg C)
484  , fbot & ! heat flux to ice bottom (W/m^2)
485  , rside ! fraction of ice that melts laterally
486 !
487 !EOP
488 !
489  integer (kind=int_kind) :: &
490  i, j & ! horizontal indices
491  , ni & ! thickness category index
492  , k & ! layer index
493  , ij & ! horizontal index, combines i and j loops
494  , icells ! number of cells with ice melting
495 
496  integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)) :: &
497  indxi, indxj ! compressed indices for cells with ice melting
498 
499  real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi) :: &
500  etot & ! total energy in column
501  , fside ! lateral heat flux (W/m^2)
502 
503  real (kind=dbl_kind) :: &
504  deltat &! SST - Tbot >= 0
505  , ustar &! skin friction velocity for fbot (m/s)
506  , wlat &! lateral melt rate (m/s)
507  , xtmp ! temporarty variable
508 
509  ! Parameters for bottom melting
510 
511  ! 0.006 = unitless param for basal heat flx ala McPhee and Maykut
512 
513  real (kind=dbl_kind), parameter :: &
514  cpchr = -cp_ocn*rhow*0.006_dbl_kind &
515  , ustar_min = 5.e-3_dbl_kind
516 
517  ! Parameters for lateral melting
518 
519  real (kind=dbl_kind), parameter :: &
520  floediam = 300.0_dbl_kind & ! effective floe diameter (m)
521  , alpha = 0.66_dbl_kind & ! constant from Steele (unitless)
522  , m1 = 1.6e-6_dbl_kind & ! constant from Maykut & Perovich
523  ! (m/s/deg^(-m2))
524  , m2 = 1.36_dbl_kind ! constant from Maykut & Perovich
525  ! (unitless)
526 
527  !-----------------------------------------------------------------
528  ! Set bottom ice surface temperature and initialize fluxes.
529  !-----------------------------------------------------------------
530 
531 
532  do j = jlo, jhi
533  do i = ilo, ihi
534  tbot(i,j) = tf(i,j)
535  fbot(i,j) = c0i
536  rside(i,j) = c0i
537  fside(i,j) = c0i
538  etot(i,j) = c0i
539  enddo ! i
540  enddo ! j
541 
542  !-----------------------------------------------------------------
543  ! Identify grid cells where ice can melt.
544  !-----------------------------------------------------------------
545 
546  icells = 0
547  do j = jlo, jhi
548  do i = ilo, ihi
549  if (aice(i,j) > puny .and. frzmlt(i,j) < c0i) then ! ice can melt
550  icells = icells + 1
551  indxi(icells) = i
552  indxj(icells) = j
553  endif
554  enddo ! i
555  enddo ! j
556 
557  do ij = 1, icells ! cells where ice can melt
558  i = indxi(ij)
559  j = indxj(ij)
560 
561  !-----------------------------------------------------------------
562  ! Use boundary layer theory for fbot.
563  ! See Maykut and McPhee (1995): JGR, 100, 24,691-24,703.
564  !-----------------------------------------------------------------
565 
566  deltat = max((sst(i,j)-tbot(i,j)),c0i)
567 
568  ! strocnx has units N/m^2 so strocnx/rho has units m^2/s^2
569  ustar = sqrt(sqrt(strocnxt(i,j)**2+strocnyt(i,j)**2)/rhow)
570  ustar = max(ustar,ustar_min*ustar_scale)
571 
572  fbot(i,j) = cpchr * deltat * ustar ! < 0
573 
574  fbot(i,j) = max(fbot(i,j), frzmlt(i,j)) ! frzmlt < fbot < 0
575 
576 !!! uncomment to use all frzmlt for standalone runs
577 ! fbot(i,j) = min (c0i, frzmlt(i,j))
578 
579  !-----------------------------------------------------------------
580  ! Compute rside. See these references:
581  ! Maykut and Perovich (1987): JGR, 92, 7032-7044
582  ! Steele (1992): JGR, 97, 17,729-17,738
583  !-----------------------------------------------------------------
584 
585  wlat = m1 * deltat**m2 ! Maykut & Perovich
586 ! rside(i,j) = wlat*dt*pi/(alpha*floediam) ! Steele
587  rside(i,j) = wlat*dtice*pi/(alpha*floediam) ! Steele
588 
589 
590  rside(i,j) = max(c0i,min(rside(i,j),c1i))
591 
592  enddo ! ij
593 
594  !-----------------------------------------------------------------
595  ! Compute heat flux associated with this value of rside.
596  !-----------------------------------------------------------------
597 
598  do ni = 1, ncat
599 
600  ! melting energy/unit area in each column, etot < 0
601  do ij = 1, icells
602  i = indxi(ij)
603  j = indxj(ij)
604  etot(i,j) = esnon(i,j,ni)
605  enddo ! ij
606 
607  do k = 1, nilyr
608 !DIR$ CONCURRENT !Cray
609 !cdir nodep !NEC
610 !ocl novrec !Fujitsu
611  do ij = 1, icells
612  i = indxi(ij)
613  j = indxj(ij)
614  etot(i,j) = etot(i,j) + eicen(i,j,ilyr1(ni)+k-1)
615  enddo ! ij
616  enddo ! k
617 
618 !DIR$ CONCURRENT !Cray
619 !cdir nodep !NEC
620 !ocl novrec !Fujitsu
621  do ij = 1, icells
622  i = indxi(ij)
623  j = indxj(ij)
624  ! lateral heat flux
625 ! fside(i,j) = fside(i,j) + rside(i,j)*etot(i,j)/dt ! fside < 0
626  fside(i,j) = fside(i,j) + rside(i,j)*etot(i,j)/dtice ! fside < 0
627  enddo ! ij
628 
629  enddo ! n
630 
631  !-----------------------------------------------------------------
632  ! Limit bottom and lateral heat fluxes if necessary.
633  !-----------------------------------------------------------------
634 
635 !DIR$ CONCURRENT !Cray
636 !cdir nodep !NEC
637 !ocl novrec !Fujitsu
638 
639 
640  do ij = 1, icells
641  i = indxi(ij)
642  j = indxj(ij)
643 
644  xtmp = frzmlt(i,j)/(fbot(i,j) + fside(i,j) + puny)
645  xtmp = min(xtmp, c1i)
646  fbot(i,j) = fbot(i,j) * xtmp
647  rside(i,j) = rside(i,j) * xtmp
648  enddo ! ij
649 
650  end subroutine frzmlt_bottom_lateral
651 
652 !=======================================================================
653 !BOP
654 !
655 ! !ROUTINE: init_thermo_vars - initialize thermodynamic variables
656 !
657 ! !DESCRIPTION:
658 !
659 ! For current thickness category, initialize the thermodynamic
660 ! variables that are aggregated and sent to the coupler, along
661 ! with other fluxes passed among subroutines.
662 !
663 ! !REVISION HISTORY:
664 !
665 ! author William H. Lipscomb, LANL
666 ! C. M. Bitz, UW
667 !
668 ! !INTERFACE:
669 !
670  subroutine init_thermo_vars(strxn, stryn, Trefn, Qrefn, &
671  & fsensn, flatn, fswabsn, flwoutn, &
672  & evapn, freshn, fsaltn, fhnetn, &
673  & fswthrun, fsurf, fcondtop, fcondbot, &
674  & fswint, einit, efinal, mvap)
675 !
676 ! !USES:
677 !
678 ! !INPUT/OUTPUT PARAMETERS:
679 !
680  real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out) :: &
681 
682  ! fields aggregated for coupler
683  strxn & ! air/ice zonal strss (N/m^2)
684  , stryn & ! air/ice merdnl strss (N/m^2)
685  , trefn & ! air tmp rfrnc level (K)
686  , qrefn & ! air sp hum rfrnc level (kg/kg)
687  , fsensn & ! surface downward sensible heat (W m-2)
688  , flatn & ! surface downward latent heat (W m-2)
689  , fswabsn & ! SW absorbed by ice (W m-2)
690  , flwoutn & ! upward LW at surface (W m-2)
691  , evapn & ! flux of vapor, atmos to ice (kg m-2 s-1)
692  , freshn & ! flux of water, ice to ocean (kg m-2 s-1)
693  , fsaltn & ! flux of salt, ice to ocean (kg m-2 s-1)
694  , fhnetn & ! fbot, corrected for surplus energy (W m-2)
695  , fswthrun & ! SW through ice to ocean (W m-2)
696 
697  ! other fields passed among subroutines
698  , fsurf & ! net flux to top surface, not including fcondtop
699  , fcondtop & ! downward cond flux at top surface (W m-2)
700  , fcondbot & ! downward cond flux at bottom surface (W m-2)
701  , fswint & ! SW absorbed in ice interior, below surface (W m-2)
702  , einit & ! initial energy of melting (J m-2)
703  , efinal & ! final energy of melting (J m-2)
704  , mvap ! ice/snow mass sublimated/condensed (kg m-2)
705 !
706 !EOP
707 !
708  integer (kind=int_kind) :: &
709  i, j ! horizontal indices
710 
711  do j = jlo, jhi
712  do i = ilo, ihi
713 
714  strxn(i,j) = c0i
715  stryn(i,j) = c0i
716  trefn(i,j) = c0i
717  qrefn(i,j) = c0i
718  fsensn(i,j) = c0i
719  flatn(i,j) = c0i
720  fswabsn(i,j) = c0i
721  flwoutn(i,j) = c0i
722  evapn(i,j) = c0i
723  freshn(i,j) = c0i
724  fsaltn(i,j) = c0i
725  fhnetn(i,j) = c0i
726  fswthrun(i,j) = c0i
727 
728  fsurf(i,j) = c0i
729  fcondtop(i,j) = c0i
730  fcondbot(i,j) = c0i
731  fswint(i,j) = c0i
732  einit(i,j) = c0i
733  efinal(i,j) = c0i
734  mvap(i,j) = c0i
735 
736  enddo ! i
737  enddo ! j
738 
739  end subroutine init_thermo_vars
740 
741 !=======================================================================
742 !BOP
743 !
744 ! !ROUTINE: init_vertical_profile - initial thickness, enthalpy, temperature
745 !
746 ! !DESCRIPTION:
747 !
748 ! Given the state variables (vicen, vsnon, eicen, esnon, Tsfcn),
749 ! compute variables needed for the vertical thermodynamics
750 ! (hin, hsn, qin, qsn, Tin, Tsn, Tsf).
751 !
752 ! !REVISION HISTORY:
753 !
754 ! authors William H. Lipscomb, LANL
755 ! C. M. Bitz, UW
756 !
757 ! !INTERFACE:
758 !
759  subroutine init_vertical_profile &
760  (ni, icells, indxi, indxj, &
761  hin, hlyr, hsn, hin_init, hsn_init, &
762  qin, tin, qsn, tsn, tsf, einit)
763 !
764 ! !USES:
765 !
766 ! use ice_exit
767 !
768 ! !INPUT/OUTPUT PARAMETERS:
769 !
770  integer (kind=int_kind), intent(in) :: &
771  ni & ! thickness category index
772  , icells ! number of cells with aicen > puny
773 
774  integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), &
775  intent(in) :: &
776  indxi, indxj ! compressed indices for cells with aicen > puny
777 
778  real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out) :: &
779  hin &! ice thickness (m)
780  , hsn &! snow thickness (m)
781  , hlyr &! ice layer thickness
782  , hin_init &! initial value of hin
783  , hsn_init &! initial value of hsn
784  , qsn &! snow enthalpy
785  , tsn &! snow temperature
786  , tsf ! ice/snow surface temperature, Tsfcn
787 
788  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nilyr), &
789  intent(out) :: &
790  qin &! ice layer enthalpy (J m-3)
791  , tin ! internal ice layer temperatures
792 
793  real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(inout) :: &
794  einit ! initial energy of melting (J m-2)
795 !
796 !EOP
797 !
798  real (kind=dbl_kind), parameter :: &
799  tmin = -100._dbl_kind ! min allowed internal temperature (deg C)
800 
801  integer (kind=int_kind) :: &
802  i, j & ! horizontal indices
803  , ij & ! horizontal index, combines i and j loops
804  , k ! ice layer index
805 
806  real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi) :: &
807  tmax ! maximum allowed snow temperature (deg C)
808 
809  real (kind=dbl_kind) :: &
810  aa1, bb1, cc1 ! terms in quadratic formula
811 
812  logical (kind=log_kind) :: & ! for vector-friendly error checks
813  tsno_high & ! flag for Tsn > Tmax
814  , tice_high & ! flag for Tin > Tmlt
815  , tsno_low & ! flag for Tsn < Tmin
816  , tice_low ! flag for Tin < Tmin
817 
818  !-----------------------------------------------------------------
819  ! Initialize arrays
820  !-----------------------------------------------------------------
821 
822  do j = jlo, jhi
823  do i = ilo, ihi
824  hin(i,j) = c0i
825  hsn(i,j) = c0i
826  hlyr(i,j) = c0i
827  hin_init(i,j) = c0i
828  hsn_init(i,j) = c0i
829  qsn(i,j) = c0i
830  tsn(i,j) = c0i
831  tsf(i,j) = c0i
832  tmax(i,j) = c0i
833  enddo
834  enddo
835 
836  do k = 1, nilyr
837  do j = jlo, jhi
838  do i = ilo, ihi
839  qin(i,j,k) = c0i
840  tin(i,j,k) = c0i
841  enddo
842  enddo
843  enddo
844 
845  !-----------------------------------------------------------------
846  ! Initialize flags
847  !-----------------------------------------------------------------
848 
849  tsno_high = .false.
850  tice_high = .false.
851  tsno_low = .false.
852  tice_low = .false.
853 
854  !-----------------------------------------------------------------
855  ! Load arrays for vertical thermo calculation.
856  !-----------------------------------------------------------------
857 !DIR$ CONCURRENT !Cray
858 !cdir nodep !NEC
859 !ocl novrec !Fujitsu
860  do ij = 1, icells
861  i = indxi(ij)
862  j = indxj(ij)
863 
864  !-----------------------------------------------------------------
865  ! Surface temperature, ice and snow thickness, snow enthalpy
866  !
867  ! Tmax based on the idea that dT ~ dq / (rhos*cp_ice)
868  ! dq ~ q dv / v
869  ! dv ~ eps11
870  ! where 'd' denotes an error due to roundoff.
871  !-----------------------------------------------------------------
872 
873  tsf(i,j) = tsfcn(i,j,ni)
874  hin(i,j) = vicen(i,j,ni) / aicen(i,j,ni)
875  hlyr(i,j) = hin(i,j) / real(nilyr,kind=dbl_kind)
876  hsn(i,j) = vsnon(i,j,ni) / aicen(i,j,ni)
877  hin_init(i,j) = hin(i,j)
878  hsn_init(i,j) = hsn(i,j)
879 
880  if (hsn(i,j) > hsnomin) then
881  qsn(i,j) = esnon(i,j,ni) / vsnon(i,j,ni) ! qsn, esnon < 0
882  tmax(i,j) = -qsn(i,j)*eps11 / (rhos*cp_ice*vsnon(i,j,ni))
883  else
884  qsn(i,j) = -rhos * lfresh
885  tmax(i,j) = puny
886  endif
887 
888  !-----------------------------------------------------------------
889  ! Compute snow temperatures from enthalpies.
890  !-----------------------------------------------------------------
891  tsn(i,j) = (lfresh + qsn(i,j)/rhos)/cp_ice ! qsn <= -rhos*Lfresh
892 
893  !-----------------------------------------------------------------
894  ! Check for Tsn > Tmax (allowing for roundoff error)
895  ! and Tsn < Tmin.
896  !-----------------------------------------------------------------
897  if (tsn(i,j) > tmax(i,j)) then
898  tsno_high = .true.
899  elseif (tsn(i,j) < tmin) then
900  tsno_low = .true.
901  endif
902 
903  enddo
904 
905  !-----------------------------------------------------------------
906  ! If Tsn is out of bounds, print diagnostics and exit.
907  !-----------------------------------------------------------------
908  if (tsno_high) then
909  do ij = 1, icells
910  i = indxi(ij)
911  j = indxj(ij)
912 
913  if (tsn(i,j) > tmax(i,j)) then ! allowing for roundoff error
914 ! write(nu_diag,*) ' '
915 ! write(nu_diag,*) 'Starting thermo, Tsn > Tmax, cat',ni
916 ! write(nu_diag,*) 'Tsn=',Tsn(i,j)
917 ! write(nu_diag,*) 'Tmax=',Tmax(i,j)
918 ! write(nu_diag,*) 'istep1, my_task, i, j:', &
919 ! istep1, my_task, i, j
920 ! write(nu_diag,*) 'qsn',qsn(i,j)
921 ! stoplabel = 'ice state at stop'
922 ! call print_state(stoplabel,i,j)
923 ! call abort_ice('vertical thermo: Tsn must be <= 0')
924  endif
925 
926  enddo ! ij
927  endif ! tsno_high
928 
929  if (tsno_low) then
930  do ij = 1, icells
931  i = indxi(ij)
932  j = indxj(ij)
933 
934 ! if (Tsn(i,j) < Tmin) then ! allowing for roundoff error
935 ! write(nu_diag,*) ' '
936 ! write(nu_diag,*) 'Starting thermo, Tsn < Tmin, cat',ni
937 ! write(nu_diag,*) 'Tsn=', Tsn(i,j)
938 ! write(nu_diag,*) 'Tmin=', Tmin
939 ! write(nu_diag,*) 'istep1, my_task, i, j:', &
940 ! istep1, my_task, i, j
941 ! write(nu_diag,*) 'qsn', qsn(i,j)
942 ! stoplabel = 'ice state at stop'
943 ! call print_state(stoplabel,i,j)
944 ! call abort_ice('vertical thermo: Tsn must be <= 0')
945 ! endif
946 
947  enddo ! ij
948  endif ! tsno_low
949 
950  !-----------------------------------------------------------------
951  ! initial energy per unit area of ice/snow, relative to 0 C
952  ! incremented later for ice
953  !-----------------------------------------------------------------
954 
955 !DIR$ CONCURRENT !Cray
956 !cdir nodep !NEC
957 !ocl novrec !Fujitsu
958  do ij = 1, icells
959  i = indxi(ij)
960  j = indxj(ij)
961 
962  if (tsn(i,j) > c0i) then ! correct roundoff error
963  tsn(i,j) = c0i
964  qsn(i,j) = -rhos*lfresh
965  endif
966 
967  einit(i,j) = hsn(i,j)*qsn(i,j)
968 
969  enddo ! ij
970 
971  do k = 1, nilyr
972 !DIR$ CONCURRENT !Cray
973 !cdir nodep !NEC
974 !ocl novrec !Fujitsu
975  do ij = 1, icells
976  i = indxi(ij)
977  j = indxj(ij)
978 
979  !-----------------------------------------------------------------
980  ! Compute ice enthalpy
981  !-----------------------------------------------------------------
982 
983  qin(i,j,k) = eicen(i,j,ilyr1(ni)+k-1) & ! qin < 0
984  * real(nilyr,kind=dbl_kind)/vicen(i,j,ni)
985 
986  !-----------------------------------------------------------------
987  ! Compute ice temperatures from enthalpies using quadratic formula
988  !-----------------------------------------------------------------
989 
990  if (l_brine) then
991  aa1 = cp_ice
992  bb1 = (cp_ocn-cp_ice)*tmlt(k) - qin(i,j,k)/rhoi - lfresh
993  cc1 = lfresh * tmlt(k)
994  tin(i,j,k) = (-bb1 - sqrt(bb1*bb1 - c4i*aa1*cc1)) / &
995  (c2i*aa1)
996  tmax(i,j) = tmlt(k)
997 
998  else ! fresh ice
999  tin(i,j,k) = (lfresh + qin(i,j,k)/rhoi) / cp_ice
1000  tmax(i,j) = -qin(i,j,k)*eps11/(rhos*cp_ice*vicen(i,j,ni))
1001  ! as above for snow
1002  endif
1003 
1004  IF (tin(i,j,k) < tmin)THEN
1005  !(Tin=Tmlt(k))
1006  eicen(i,j,ilyr1(ni)+k-1) = &
1007  (rhoi *cp_ocn*tmlt(k)) &
1008  * vicen(i,j,ni)/real(nilyr,kind=dbl_kind)
1009  qin(i,j,k) = eicen(i,j,ilyr1(ni)+k-1) & ! qin < 0
1010  * real(nilyr,kind=dbl_kind)/vicen(i,j,ni)
1011  aa1 = cp_ice
1012  bb1 = (cp_ocn-cp_ice)*tmlt(k) - qin(i,j,k)/rhoi - lfresh
1013  cc1 = lfresh * tmlt(k)
1014  tin(i,j,k) = (-bb1 - sqrt(bb1*bb1 - c4i*aa1*cc1)) / &
1015  (c2i*aa1)
1016  end if
1017 
1018 
1019 
1020  !-----------------------------------------------------------------
1021  ! Check for Tin > Tmax and Tin < Tmin
1022  !-----------------------------------------------------------------
1023 
1024  if (tin(i,j,k) > tmax(i,j)) then
1025  tice_high = .true.
1026  elseif (tin(i,j,k) < tmin) then
1027  tice_low = .true.
1028  endif
1029 
1030  enddo ! ij
1031 
1032  !-----------------------------------------------------------------
1033  ! If Tin is out of bounds, print diagnostics and exit.
1034  !-----------------------------------------------------------------
1035 
1036  if (tice_high) then
1037  do ij = 1, icells
1038  i = indxi(ij)
1039  j = indxj(ij)
1040 
1041  if (tin(i,j,k) > tmax(i,j)) then
1042 ! write(nu_diag,*) ' '
1043 ! write(nu_diag,*) 'Starting thermo, T > Tmax, cat',&
1044 ! ni,', layer',k
1045 ! write(nu_diag,*) 'Tin=',Tin(i,j,k),', Tmax=',Tmax(i,j)
1046 ! write(nu_diag,*) 'istep1, my_task, i, j:', &
1047 ! istep1, my_task, i, j
1048 ! write(nu_diag,*) 'qin',qin(i,j,k)
1049 ! stoplabel = 'ice state at stop'
1050 ! call print_state(stoplabel,i,j)
1051 ! call abort_ice('vertical thermo: Tin must be <= Tmax')
1052  endif
1053 
1054  enddo ! ij
1055  endif ! tice_high
1056 
1057  if (tice_low) then
1058  do ij = 1, icells
1059  i = indxi(ij)
1060  j = indxj(ij)
1061 
1062 ! if (Tin(i,j,k) < Tmin) then
1063 ! write(nu_diag,*) ' '
1064 ! write(nu_diag,*) 'Starting thermo T < Tmin, cat', &
1065 ! ni,', layer',k
1066 ! write(nu_diag,*) 'Tin =', Tin(i,j,k)
1067 ! write(nu_diag,*) 'Tmin =', Tmin
1068 !! write(nu_diag,*) 'istep1, my_task, i, j:', &
1069 ! istep1, my_task, i, j
1070 ! stoplabel = 'ice state at stop'
1071 ! call print_state(stoplabel,i,j)
1072 ! call abort_ice('vertical_thermo: Tin must be > Tmin')
1073 ! endif
1074  enddo ! ij
1075  endif ! tice_low
1076 
1077  !-----------------------------------------------------------------
1078  ! initial energy per unit area of ice/snow, relative to 0 C
1079  !-----------------------------------------------------------------
1080 
1081 !DIR$ CONCURRENT !Cray
1082 !cdir nodep !NEC
1083 !ocl novrec !Fujitsu
1084  do ij = 1, icells
1085  i = indxi(ij)
1086  j = indxj(ij)
1087 
1088  if (tin(i,j,k) > c0i) then ! correct roundoff error
1089  tin(i,j,k) = c0i
1090  qin(i,j,k) = -rhoi*lfresh
1091  endif
1092 
1093  einit(i,j) = einit(i,j) + hlyr(i,j)*qin(i,j,k)
1094 
1095  enddo ! ij
1096  enddo ! k
1097 
1098  end subroutine init_vertical_profile
1099 
1100 !=======================================================================
1101 !BOP
1102 !
1103 ! !ROUTINE: temperature_changes - new vertical temperature profile
1104 !
1105 ! !DESCRIPTION:
1106 !
1107 ! Compute new surface temperature and internal ice and snow
1108 ! temperatures. Include effects of salinity on sea ice heat
1109 ! capacity in a way that conserves energy (Bitz and Lipscomb, 1999).
1110 !
1111 ! New temperatures are computed iteratively by solving a tridiagonal
1112 ! system of equations; heat capacity is updated with each iteration.
1113 ! Finite differencing is backward implicit.
1114 !
1115 ! !REVISION HISTORY:
1116 !
1117 ! authors William H. Lipscomb, LANL
1118 ! C. M. Bitz, UW
1119 !
1120 ! !INTERFACE:
1121 !
1122  subroutine temperature_changes &
1123  (ni, icells, indxi, indxj, &
1124  hlyr, hsn, qin, tin, qsn, tsn, &
1125  tsf, tbot, fbot, fsensn, flatn, fswabsn, &
1126  flwoutn, fswthrun, fhnetn, fsurf, fcondtop, fcondbot, &
1127  fswint, einit)
1128 !
1129 ! !USES:
1130 !
1131 ! use ice_exit
1132 !
1133 ! !INPUT/OUTPUT PARAMETERS:
1134 !
1135  integer (kind=int_kind), intent(in) :: &
1136  ni & ! thickness category index
1137  , icells ! number of cells with aicen > puny
1138 
1139  integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), &
1140  intent(in) :: &
1141  indxi, indxj ! compressed indices for cells with aicen > puny
1142 
1143  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout):: &
1144  hlyr & ! ice layer thickness
1145  , hsn & ! snow thickness (m)
1146  , tbot & ! ice bottom surface temperature (deg C)
1147  , fbot & ! ice-ocean heat flux at bottom surface (W/m^2)
1148  , qsn & ! snow enthalpy
1149  , tsn & ! internal snow temperature
1150  , tsf & ! ice/snow surface temperature, Tsfcn
1151  , fsensn & ! surface downward sensible heat (W m-2)
1152  , flatn & ! surface downward latent heat (W m-2)
1153  , fswabsn & ! shortwave absorbed by ice (W m-2)
1154  , flwoutn & ! upward LW at surface (W m-2)
1155  , fswthrun & ! SW through ice to ocean (W m-2)
1156  , fhnetn & ! fbot, corrected for any surplus energy
1157  , fsurf & ! net flux to top surface, not including fcondtop
1158  , fcondtop & ! downward cond flux at top surface (W m-2)
1159  , fcondbot & ! downward cond flux at bottom surface (W m-2)
1160  , fswint ! SW absorbed in ice interior, below surface (W m-2)
1161 
1162  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nilyr), &
1163  intent(inout) :: &
1164  qin & ! ice layer enthalpy (J m-3)
1165  , tin ! internal ice layer temperatures
1166 
1167  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(in):: &
1168  einit ! initial energy of melting (J m-2)
1169 !
1170 !EOP
1171 !
1172  integer (kind=int_kind), parameter :: &
1173  nitermax = 50 ! max number of iterations in temperature solver
1174 
1175  real (kind=dbl_kind), parameter :: &
1176  alph = c3i & ! constant used to get 2nd order accurate fluxes
1177  , bet = -p333 & ! constant used to get 2nd order accurate fluxes
1178  , tsf_errmax = 5.e-4 ! max allowed error in Tsf
1179  ! recommend Tsf_errmax < 0.01 K
1180 
1181  integer (kind=int_kind) :: &
1182  i, j & ! horizontal indices
1183  , ij & ! horizontal index, combines i and j loops
1184  , k & ! ice layer index
1185  , nmat & ! matrix dimension
1186  , niter & ! iteration counter in temperature solver
1187  , isolve ! number of cells with temps not converged
1188 
1189  integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)) :: &
1190  indxii, indxjj ! compressed indices for cells not converged
1191 
1192  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi) :: &
1193  tsn_init & ! Tsn at beginning of time step
1194  , tsn_start & ! Tsn at start of iteration
1195  , tsf_start & ! Tsf at start of iteration
1196  , dtsf & ! Tsf - Tsf_start
1197  , dtsf_prev & ! dTsf from previous iteration
1198  , dfsurf_dt & ! derivative of fsurf wrt Tsf
1199  , dfsens_dt & ! deriv of fsens wrt Tsf (W m-2 deg-1)
1200  , dflat_dt & ! deriv of flat wrt Tsf (W m-2 deg-1)
1201  , dflwout_dt & ! deriv of flwout wrt Tsf (W m-2 deg-1)
1202  , fswsfc ! SW absorbed at ice/snow surface (W m-2)
1203 
1204  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi) :: &
1205  khs & ! ksno / hsn
1206  , etas & ! dt / (rhos * cp_ice * hsn)
1207  , dt_rhoi_hlyr & ! dt/(rhoi*hlyr)
1208  , avg_tin & ! = 1. if Tin averaged w/Tin_start, else = 0.
1209  , enew ! new energy of melting after temp change (J m-2)
1210 
1211  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nilyr) :: &
1212  tin_init & ! Tin at beginning of time step
1213  , tin_start & ! Tin at start of iteration
1214  , iabs & ! SW absorbed in particular layer (W m-2)
1215  , etai ! dt / (rho * cp * h) for ice layers
1216 
1217  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nilyr+1) :: &
1218  khi ! ki / hlyr
1219 
1220  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nilyr+2) :: &
1221  diag & ! diagonal matrix elements
1222  , sbdiag & ! sub-diagonal matrix elements
1223  , spdiag & ! super-diagonal matrix elements
1224  , rhs & ! rhs of tri-diagonal matrix equation
1225  , tmat ! matrix output temperatures
1226 
1227  real (kind=dbl_kind) :: &
1228  ci & ! specific heat of sea ice (J kg-1 deg-1)
1229  , spdiag2 & ! term to right of superdiag on top row
1230  , avg_tsf & ! = 1. if Tsf averaged w/Tsf_start, else = 0.
1231  , avg_tsn & ! = 1. if Tsn averaged w/Tsn_start, else = 0.
1232  , ferr & ! energy conservation error (W m-2)
1233  , w1,w2,w3,w4,w5,w6,w7 ! temporary variables
1234 
1235  logical (kind=log_kind), dimension (ilo:ihi,jlo:jhi) :: &
1236  converged ! = true when local solution has converged
1237 
1238  logical (kind=log_kind) :: &
1239  all_converged ! = true when all cells have converged
1240 
1241  !-----------------------------------------------------------------
1242  ! Initialize
1243  !-----------------------------------------------------------------
1244 
1245  all_converged = .false.
1246 
1247 !DIR$ CONCURRENT !Cray
1248 !cdir nodep !NEC
1249 !ocl novrec !Fujitsu
1250  do ij = 1, icells
1251 
1252  i = indxi(ij)
1253  j = indxj(ij)
1254 
1255  converged(i,j) = .false.
1256  khs(i,j) = c0i
1257  dtsf(i,j) = c0i
1258  dtsf_prev(i,j) = c0i
1259  tsf_start(i,j) = c0i
1260  enew(i,j) = c0i
1261 
1262  dfsurf_dt(i,j) = c0i
1263  dfsens_dt(i,j) = c0i
1264  dflat_dt(i,j) = c0i
1265  dflwout_dt(i,j) = c0i
1266 
1267  fhnetn(i,j) = fbot(i,j) ! ocean energy used by the ice, <= 0
1268 
1269  fswsfc(i,j) = c0i
1270 ! dt_rhoi_hlyr(i,j) = dt / (rhoi*hlyr(i,j)) ! used by matrix solver
1271  dt_rhoi_hlyr(i,j) = dtice / (rhoi*hlyr(i,j)) ! used by matrix solver
1272 
1273  if (hsn(i,j) > hsnomin) then
1274 ! etas(i,j) = dt / (rhos*cp_ice*hsn(i,j)) ! used by matrix solver
1275  etas(i,j) = dtice / (rhos*cp_ice*hsn(i,j)) ! used by matrix solver
1276  else
1277  etas(i,j) = c0i
1278  endif
1279 
1280  tsn_init(i,j) = tsn(i,j) ! beginning of time step
1281  tsn_start(i,j) = tsn(i,j) ! beginning of iteration
1282 
1283  enddo ! ij
1284 
1285  do k = 1, nilyr
1286 !DIR$ CONCURRENT !Cray
1287 !cdir nodep !NEC
1288 !ocl novrec !Fujitsu
1289  do ij = 1, icells
1290  i = indxi(ij)
1291  j = indxj(ij)
1292 
1293  iabs(i,j,k) = c0i
1294  tin_init(i,j,k) = tin(i,j,k) ! beginning of time step
1295  tin_start(i,j,k) = tin(i,j,k) ! beginning of iteration
1296  enddo ! ij
1297  enddo ! k
1298 
1299  !-----------------------------------------------------------------
1300  ! Compute thermal conductivity at interfaces (held fixed during
1301  ! the subsequent iteration).
1302  !-----------------------------------------------------------------
1303  call conductivity &
1304  (icells, indxi, indxj, &
1305  hlyr, hsn, tin, tbot, khi, khs)
1306 
1307  !-----------------------------------------------------------------
1308  ! Compute solar radiation absorbed in ice and penetrating to ocean
1309  !-----------------------------------------------------------------
1310  call absorbed_solar &
1311  (ni, icells, indxi, indxj, &
1312  hlyr, hsn, fswsfc, fswint, fswthrun, iabs)
1313 
1314  !-----------------------------------------------------------------
1315  ! Solve for new temperatures.
1316  ! Iterate until temperatures converge with minimal energy error.
1317  !-----------------------------------------------------------------
1318 
1319  do niter = 1, nitermax
1320 
1321  if (all_converged) then ! thermo calculation is done
1322  exit
1323  else ! identify cells not yet converged
1324  isolve = 0
1325  do ij = 1, icells
1326  i = indxi(ij)
1327  j = indxj(ij)
1328  if (.not.converged(i,j)) then
1329  isolve = isolve + 1
1330  indxii(isolve) = i
1331  indxjj(isolve) = j
1332  endif
1333  enddo ! ij
1334  endif
1335 
1336  !-----------------------------------------------------------------
1337  ! Update radiative and turbulent fluxes and their derivatives
1338  ! with respect to Tsf.
1339  !-----------------------------------------------------------------
1340 
1341  call surface_fluxes &
1342  (isolve, indxii, indxjj, &
1343  tsf, fswsfc, &
1344  flwoutn, fsensn, flatn, fsurf, &
1345  dflwout_dt, dfsens_dt, dflat_dt, dfsurf_dt)
1346 
1347 !DIR$ CONCURRENT !Cray
1348 !cdir nodep !NEC
1349 !ocl novrec !Fujitsu
1350  do ij = 1, isolve
1351  i = indxii(ij) ! NOTE: not indxi and indxj
1352  j = indxjj(ij)
1353 
1354  !-----------------------------------------------------------------
1355  ! Compute conductive flux at top surface, fcondtop.
1356  ! If fsurf < fcondtop and Tsf = 0, then reset Tsf to slightly less
1357  ! than zero (but not less than -puny).
1358  !-----------------------------------------------------------------
1359 
1360  if (hsn(i,j) > hsnomin) then
1361  fcondtop(i,j) = c2i * khs(i,j) * (tsf(i,j) - tsn(i,j))
1362  else
1363  fcondtop(i,j) = khi(i,j,1) &
1364  * (alph*(tsf(i,j) - tin(i,j,1)) &
1365  + bet*(tsf(i,j) - tin(i,j,2)))
1366  endif
1367  if (fsurf(i,j) < fcondtop(i,j)) &
1368  tsf(i,j) = min(tsf(i,j), -puny)
1369 
1370  !-----------------------------------------------------------------
1371  ! Save surface temperature at start of iteration
1372  !-----------------------------------------------------------------
1373  tsf_start(i,j) = tsf(i,j)
1374 
1375  enddo ! ij
1376 
1377  !-----------------------------------------------------------------
1378  ! Compute specific heat of each ice layer
1379  !-----------------------------------------------------------------
1380 
1381  do k = 1, nilyr
1382  do ij = 1, isolve
1383  i = indxii(ij)
1384  j = indxjj(ij)
1385 
1386  if (l_brine) then
1387  ci = cp_ice - lfresh*tmlt(k) / &
1388  (tin_init(i,j,k)*tin(i,j,k))
1389  else
1390  ci = cp_ice
1391  endif
1392  etai(i,j,k) = dt_rhoi_hlyr(i,j) / ci
1393 
1394  enddo ! ij
1395  enddo ! k
1396 
1397 !DIR$ CONCURRENT !Cray
1398 !cdir nodep !NEC
1399 !ocl novrec !Fujitsu
1400  do ij = 1, isolve
1401  i = indxii(ij)
1402  j = indxjj(ij)
1403 
1404  !-----------------------------------------------------------------
1405  ! Compute matrix elements
1406  !
1407  ! Four possible cases to solve:
1408  ! (1) Cold surface (Tsf < 0), snow present
1409  ! (2) Melting surface (Tsf = 0), snow present
1410  ! (3) Cold surface (Tsf < 0), no snow
1411  ! (4) Melting surface (Tsf = 0), no snow
1412  !-----------------------------------------------------------------
1413 
1414  !-----------------------------------------------------------------
1415  ! first row
1416  ! Tsf for case 1; dummy equation for cases 2, 3 and 4
1417  !-----------------------------------------------------------------
1418 
1419  sbdiag(i,j,1) = c0i
1420  spdiag2 = c0i
1421 
1422  if (hsn(i,j) > hsnomin) then
1423  if (tsf(i,j) <= -puny) then
1424  spdiag(i,j,1) = c2i*khs(i,j)
1425  diag(i,j,1) = dfsurf_dt(i,j) - c2i*khs(i,j)
1426  rhs(i,j,1) = dfsurf_dt(i,j)*tsf(i,j) - fsurf(i,j)
1427  else
1428  spdiag(i,j,1) = c0i
1429  diag(i,j,1) = c1i
1430  rhs(i,j,1) = c0i
1431  endif
1432  else ! hsn <= hsnomin
1433  if (tsf(i,j) <= -puny) then
1434  spdiag(i,j,1) = c0i
1435  diag(i,j,1) = c1i
1436  rhs(i,j,1) = c0i
1437  else
1438  spdiag(i,j,1) = c0i
1439  diag(i,j,1) = c1i
1440  rhs(i,j,1) = c0i
1441  endif
1442  endif
1443 
1444  !-----------------------------------------------------------------
1445  ! second row
1446  ! Tsn for cases 1 and 2; Tsf for case 3; dummy for case 4
1447  !-----------------------------------------------------------------
1448 
1449  if (hsn(i,j) > hsnomin) then
1450  if (tsf(i,j) <= -puny) then
1451  sbdiag(i,j,2) = -etas(i,j) * c2i * khs(i,j)
1452  spdiag(i,j,2) = -etas(i,j) * khi(i,j,1)
1453  spdiag2 = c0i
1454  diag(i,j,2) = c1i &
1455  + etas(i,j) * (c2i*khs(i,j) + khi(i,j,1))
1456  rhs(i,j,2) = tsn_init(i,j)
1457  else
1458  sbdiag(i,j,2) = c0i
1459  spdiag(i,j,2) = -etas(i,j) * khi(i,j,1)
1460  spdiag2 = c0i
1461  diag(i,j,2) = c1i &
1462  + etas(i,j) * (c2i*khs(i,j) + khi(i,j,1))
1463  rhs(i,j,2) = tsn_init(i,j) &
1464  + etas(i,j)*c2i*khs(i,j)*tsf(i,j)
1465  endif
1466  else ! hsn <= hsnomin
1467  if (tsf(i,j) <= -puny) then
1468  sbdiag(i,j,2) = c0i
1469  spdiag(i,j,2) = alph * khi(i,j,1)
1470  spdiag2 = bet * khi(i,j,1)
1471  diag(i,j,2) = dfsurf_dt(i,j) - alph*khi(i,j,1) &
1472  - bet*khi(i,j,1)
1473  rhs(i,j,2) = dfsurf_dt(i,j)*tsf(i,j) - fsurf(i,j)
1474  else
1475  sbdiag(i,j,2) = c0i
1476  spdiag(i,j,2) = c0i
1477  spdiag2 = c0i
1478  diag(i,j,2) = c1i
1479  rhs(i,j,2) = c0i
1480  endif
1481  endif
1482 
1483  !-----------------------------------------------------------------
1484  ! third row, Tin(i,j,1)
1485  !
1486  ! For each internal ice layer compute the specific heat ci,
1487  ! a function of the starting temperature and of the latest
1488  ! guess for the final temperature.
1489  !-----------------------------------------------------------------
1490 
1491  if (hsn(i,j) > hsnomin) then
1492  if (tsf(i,j) <= -puny) then
1493  sbdiag(i,j,3) = -etai(i,j,1) * khi(i,j,1)
1494  spdiag(i,j,3) = -etai(i,j,1) * khi(i,j,2)
1495  diag(i,j,3) = c1i &
1496  + etai(i,j,1)*(khi(i,j,1) + khi(i,j,2))
1497  rhs(i,j,3) = tin_init(i,j,1) &
1498  + etai(i,j,1)*iabs(i,j,1)
1499  else
1500  sbdiag(i,j,3) = -etai(i,j,1) * khi(i,j,1)
1501  spdiag(i,j,3) = -etai(i,j,1) * khi(i,j,2)
1502  diag(i,j,3) = c1i &
1503  + etai(i,j,1)*(khi(i,j,1) + khi(i,j,2))
1504  rhs(i,j,3) = tin_init(i,j,1) &
1505  + etai(i,j,1)*iabs(i,j,1)
1506  endif
1507  else ! hsn <= hsnomin
1508  if (tsf(i,j) <= -puny) then
1509  sbdiag(i,j,3) = -(alph+bet) * etai(i,j,1) * khi(i,j,1)
1510  spdiag(i,j,3) = -etai(i,j,1) &
1511  * (khi(i,j,2) - bet*khi(i,j,1))
1512  diag(i,j,3) = c1i + etai(i,j,1) &
1513  * (khi(i,j,2) + alph*khi(i,j,1))
1514  rhs(i,j,3) = tin_init(i,j,1) &
1515  + etai(i,j,1)*iabs(i,j,1)
1516  else
1517  sbdiag(i,j,3) = c0i
1518  spdiag(i,j,3) = -etai(i,j,1) &
1519  * (khi(i,j,2) - bet*khi(i,j,1))
1520  diag(i,j,3) = c1i + etai(i,j,1) &
1521  * (khi(i,j,2) + alph*khi(i,j,1))
1522  rhs(i,j,3) = tin_init(i,j,1) &
1523  + etai(i,j,1)*iabs(i,j,1)
1524 !!! & + (alph+bet)*etai(i,j,k)*khi(i,j,1)*Tsf(i,j)
1525  ! commented line needed if Tsf were in Kelvin
1526  endif
1527  endif
1528 
1529  if (hsn(i,j) <= hsnomin .and. tsf(i,j) <= -puny) then ! case 3
1530  ! remove spdiag2 in row 2 by adding rows 2 and 3
1531  diag(i,j,2) = spdiag(i,j,3) * diag(i,j,2) &
1532  - spdiag2 * sbdiag(i,j,3)
1533  spdiag(i,j,2) = spdiag(i,j,3) * spdiag(i,j,2) &
1534  - spdiag2 * diag(i,j,3)
1535  rhs(i,j,2) = spdiag(i,j,3) * rhs(i,j,2) &
1536  - spdiag2 * rhs(i,j,3)
1537  endif
1538 
1539  !-----------------------------------------------------------------
1540  ! Bottom row, Tin(i,j,nilyr)
1541  !-----------------------------------------------------------------
1542 
1543  k = nilyr
1544  sbdiag(i,j,k+2) = -etai(i,j,k) &
1545  * (khi(i,j,k) - bet*khi(i,j,k+1))
1546  spdiag(i,j,k+2) = c0i
1547  diag(i,j,k+2) = c1i + etai(i,j,k) &
1548  * (khi(i,j,k) + alph*khi(i,j,k+1))
1549  rhs(i,j,k+2) = tin_init(i,j,k) + etai(i,j,k)*iabs(i,j,k) &
1550  + (alph+bet)*etai(i,j,k) &
1551  *khi(i,j,k+1)*tbot(i,j)
1552 
1553  enddo ! ij
1554 
1555  !-----------------------------------------------------------------
1556  ! Ice interior
1557  !-----------------------------------------------------------------
1558  do k = 2, nilyr-1
1559  do ij = 1, isolve
1560  i = indxii(ij)
1561  j = indxjj(ij)
1562 
1563  sbdiag(i,j,k+2) = -etai(i,j,k) * khi(i,j,k)
1564  spdiag(i,j,k+2) = -etai(i,j,k) * khi(i,j,k+1)
1565  diag(i,j,k+2) = c1i - sbdiag(i,j,k+2) - spdiag(i,j,k+2)
1566  rhs(i,j,k+2) = tin_init(i,j,k) &
1567  + etai(i,j,k) * iabs(i,j,k)
1568 
1569  enddo ! ij
1570  enddo ! k
1571 
1572  !-----------------------------------------------------------------
1573  ! Solve tridiagonal matrix to obtain the new temperatures.
1574  !-----------------------------------------------------------------
1575 
1576  nmat = nilyr + 2
1577 
1578  call tridiag_solver &
1579  (isolve, indxii, indxjj, &
1580  nmat, sbdiag, diag, spdiag, rhs, tmat)
1581 
1582  !-----------------------------------------------------------------
1583  ! Reload temperatures from matrix solution.
1584  !-----------------------------------------------------------------
1585 
1586  do ij = 1, isolve
1587  i = indxii(ij)
1588  j = indxjj(ij)
1589 
1590  if (hsn(i,j) > hsnomin) then
1591  if (tsf(i,j) <= -puny) then
1592  tsf(i,j) = tmat(i,j,1)
1593  tsn(i,j) = tmat(i,j,2)
1594  else
1595  tsf(i,j) = c0i
1596  tsn(i,j) = tmat(i,j,2)
1597  endif
1598  else ! hsn <= hsnomin
1599  if (tsf(i,j) <= -puny) then
1600  tsf(i,j) = tmat(i,j,2)
1601  else
1602  tsf(i,j) = c0i
1603  endif
1604  endif
1605 
1606  enddo
1607 
1608  do k = 1, nilyr
1609  do ij = 1, isolve
1610  i = indxii(ij)
1611  j = indxjj(ij)
1612 
1613  tin(i,j,k) = tmat(i,j,k+2)
1614  enddo
1615  enddo
1616 
1617  !-----------------------------------------------------------------
1618  ! Determine whether the computation has converged to an acceptable
1619  ! solution. Five conditions must be satisfied:
1620  !
1621  ! (1) Tsf <= 0 C.
1622  ! (2) Tsf is not oscillating; i.e., if both dTsf(niter) and
1623  ! dTsf(niter-1) have magnitudes greater than puny, then
1624  ! dTsf(niter)/dTsf(niter-1) cannot be a negative number
1625  ! with magnitude greater than 0.5.
1626  ! (3) abs(dTsf) < Tsf_errmax
1627  ! (4) If Tsf = 0 C, then the downward turbulent/radiative
1628  ! flux, fsurf, must be greater than or equal to the downward
1629  ! conductive flux, fcondtop.
1630  ! (5) The net energy added to the ice per unit time must equal
1631  ! the net change in internal ice energy per unit time,
1632  ! within the prescribed error ferrmax.
1633  !
1634  ! For briny ice (the standard case), Tsn and Tin are limited
1635  ! to prevent them from exceeding their melting temperatures.
1636  ! (Note that the specific heat formula for briny ice assumes
1637  ! that T < Tmlt.)
1638  ! For fresh ice there is no limiting, since there are cases
1639  ! when the only convergent solution has Tsn > 0 and/or Tin > 0.
1640  ! Above-zero temperatures are then reset to zero (with melting
1641  ! to conserve energy) in the thickness_changes subroutine.
1642  !-----------------------------------------------------------------
1643 
1644  ! initialize global convergence flag
1645  all_converged = .true.
1646 
1647 !DIR$ CONCURRENT !Cray
1648 !cdir nodep !NEC
1649 !ocl novrec !Fujitsu
1650  do ij = 1, isolve
1651  i = indxii(ij)
1652  j = indxjj(ij)
1653 
1654  if (l_brine) tsn(i,j) = min(tsn(i,j), c0i)
1655 
1656  !-----------------------------------------------------------------
1657  ! Initialize convergence flag (true until proven false), dTsf,
1658  ! and temperature-averaging coefficients.
1659  ! Average only if test 1 or 2 fails.
1660  !-----------------------------------------------------------------
1661 
1662  converged(i,j) = .true.
1663  dtsf(i,j) = tsf(i,j) - tsf_start(i,j)
1664  avg_tsf = c0i
1665  avg_tsn = c0i
1666  avg_tin(i,j) = c0i
1667 
1668  !-----------------------------------------------------------------
1669  ! Condition 1: check for Tsf > 0
1670  ! If Tsf > 0, set Tsf = 0, then average Tsn and Tin to force
1671  ! internal temps below their melting temps.
1672  !-----------------------------------------------------------------
1673 
1674  if (tsf(i,j) > puny) then
1675  tsf(i,j) = c0i
1676  dtsf(i,j) = -tsf_start(i,j)
1677  if (l_brine) then ! average with starting temp
1678  avg_tsn = c1i
1679  avg_tin(i,j) = c1i
1680  endif
1681  converged(i,j) = .false.
1682  all_converged = .false.
1683 
1684  !-----------------------------------------------------------------
1685  ! Condition 2: check for oscillating Tsf
1686  ! If oscillating, average all temps to increase rate of convergence.
1687  !-----------------------------------------------------------------
1688 
1689  elseif (niter > 1 &! condition (2)
1690  .and. tsf_start(i,j) <= -puny &
1691  .and. abs(dtsf(i,j)) > puny &
1692  .and. abs(dtsf_prev(i,j)) > puny &
1693  .and. -dtsf(i,j)/(dtsf_prev(i,j)+puny*puny) > p5) then
1694 
1695  if (l_brine) then ! average with starting temp
1696  avg_tsf = c1i
1697  avg_tsn = c1i
1698  avg_tin(i,j) = c1i
1699  endif
1700  dtsf(i,j) = p5 * dtsf(i,j)
1701  converged(i,j) = .false.
1702  all_converged = .false.
1703  endif
1704 
1705  !-----------------------------------------------------------------
1706  ! If condition 1 or 2 failed, average new surface/snow
1707  ! temperatures with their starting values.
1708  ! (No change if both tests passed)
1709  !-----------------------------------------------------------------
1710  tsf(i,j) = tsf(i,j) &
1711  + avg_tsf * p5 * (tsf_start(i,j) - tsf(i,j))
1712  tsn(i,j) = tsn(i,j) &
1713  + avg_tsn * p5 * (tsn_start(i,j) - tsn(i,j))
1714 
1715  !-----------------------------------------------------------------
1716  ! Compute qsn and increment new energy.
1717  !-----------------------------------------------------------------
1718  qsn(i,j) = -rhos * (lfresh - cp_ice*tsn(i,j))
1719  enew(i,j) = hsn(i,j) * qsn(i,j)
1720 
1721  tsn_start(i,j) = tsn(i,j) ! for next iteration
1722 
1723  enddo ! ij
1724 
1725  do k = 1, nilyr
1726 !DIR$ CONCURRENT !Cray
1727 !cdir nodep !NEC
1728 !ocl novrec !Fujitsu
1729  do ij = 1, isolve
1730  i = indxii(ij)
1731  j = indxjj(ij)
1732 
1733  if (l_brine) tin(i,j,k) = min(tin(i,j,k), tmlt(k))
1734 
1735  !-----------------------------------------------------------------
1736  ! If condition 1 or 2 failed, average new ice layer
1737  ! temperatures with their starting values.
1738  !-----------------------------------------------------------------
1739  tin(i,j,k) = tin(i,j,k) &
1740  + avg_tin(i,j)*p5*(tin_start(i,j,k)-tin(i,j,k))
1741 
1742  !-----------------------------------------------------------------
1743  ! Compute qin and increment new energy.
1744  !-----------------------------------------------------------------
1745  qin(i,j,k) = -rhoi * (cp_ice*(tmlt(k)-tin(i,j,k)) &
1746  + lfresh*(c1i-tmlt(k)/tin(i,j,k)) &
1747  - cp_ocn*tmlt(k))
1748  enew(i,j) = enew(i,j) + hlyr(i,j)*qin(i,j,k)
1749 
1750  tin_start(i,j,k) = tin(i,j,k) ! for next iteration
1751 
1752  enddo ! ij
1753  enddo ! k
1754 
1755 !DIR$ CONCURRENT !Cray
1756 !cdir nodep !NEC
1757 !ocl novrec !Fujitsu
1758  do ij = 1, isolve
1759  i = indxii(ij)
1760  j = indxjj(ij)
1761 
1762  !-----------------------------------------------------------------
1763  ! Condition 3: check for large change in Tsf
1764  !-----------------------------------------------------------------
1765 
1766  if (abs(dtsf(i,j)) > tsf_errmax) then
1767  converged(i,j) = .false.
1768  all_converged = .false.
1769  endif
1770 
1771  !-----------------------------------------------------------------
1772  ! Condition 4: check for fsurf < fcondtop with Tsf > 0
1773  !-----------------------------------------------------------------
1774 
1775  fsurf(i,j) = fsurf(i,j) + dtsf(i,j)*dfsurf_dt(i,j)
1776  if (hsn(i,j) > hsnomin) then
1777  fcondtop(i,j) = c2i * khs(i,j) * (tsf(i,j)-tsn(i,j))
1778  else
1779  fcondtop(i,j) = khi(i,j,1) &
1780  * (alph*(tsf(i,j)-tin(i,j,1)) &
1781  + bet*(tsf(i,j)-tin(i,j,2)))
1782  endif
1783 
1784  if (tsf(i,j) > -puny .and. fsurf(i,j) < fcondtop(i,j)) then
1785  converged(i,j) = .false.
1786  all_converged = .false.
1787  endif
1788 
1789  !-----------------------------------------------------------------
1790  ! Condition 5: check for energy conservation error
1791  ! Change in internal ice energy should equal net energy input.
1792  !-----------------------------------------------------------------
1793 
1794  fcondbot(i,j) = khi(i,j,nilyr+1) * &
1795  (alph*(tin(i,j,nilyr) - tbot(i,j)) &
1796  + bet*(tin(i,j,nilyr-1) - tbot(i,j)))
1797 
1798 ! ferr = abs( (enew(i,j)-einit(i,j))/dt &
1799  ferr = abs( (enew(i,j)-einit(i,j))/dtice &
1800  - (fcondtop(i,j) - fcondbot(i,j) + fswint(i,j)) )
1801 
1802  ! factor of 0.9 allows for roundoff errors later
1803  if (ferr > 0.9_dbl_kind*ferrmax) then ! condition (5)
1804  converged(i,j) = .false.
1805  all_converged = .false.
1806  endif
1807 
1808  dtsf_prev(i,j) = dtsf(i,j)
1809 
1810  enddo ! ij
1811 
1812  enddo ! temperature iteration
1813 
1814  if (.not.all_converged .and.dbg_set(dbg_sbrio) ) then
1815  do ij = 1, icells
1816  i = indxi(ij)
1817  j = indxj(ij)
1818 
1819  !-----------------------------------------------------------------
1820  ! Check for convergence failures.
1821  !-----------------------------------------------------------------
1822  if (.not.converged(i,j)) then
1823 ! write(nu_diag,*) 'Thermo iteration does not converge,', &
1824 ! 'istep1, my_task, i, j, n:', &
1825 ! istep1, my_task, i, j, ni
1826 ! write(nu_diag,*) 'Ice thickness:', hlyr(i,j)*nilyr
1827 ! write(nu_diag,*) 'Snow thickness:', hsn(i,j)
1828 ! write(nu_diag,*) 'dTsf, Tsf_errmax:',dTsf(i,j),Tsf_errmax
1829 ! write(nu_diag,*) 'Tsf:', Tsf(i,j)
1830 ! write(nu_diag,*) 'fsurf:', fsurf(i,j)
1831 ! write(nu_diag,*) 'fcondtop, fcondbot, fswint', &
1832 ! fcondtop(i,j), fcondbot(i,j), fswint(i,j)
1833 ! write(nu_diag,*) 'Flux conservation error =', ferr
1834 ! write(nu_diag,*) 'Initial snow and ice temperatures:'
1835 ! write(nu_diag,*) Tsn_init(i,j), &
1836 ! (Tin_init(i,j,k),k=1,nilyr)
1837 ! write(nu_diag,*) 'Final temperatures:'
1838 ! write(nu_diag,*) Tsn(i,j), (Tin(i,j,k),k=1,nilyr)
1839 ! stoplabel = 'ice state at stop'
1840 ! call print_state(stoplabel,i,j)
1841 ! call abort_ice('vertical thermo: convergence error')
1842  endif
1843  enddo ! ij
1844  endif ! all_converged
1845 
1846 !DIR$ CONCURRENT !Cray
1847 !cdir nodep !NEC
1848 !ocl novrec !Fujitsu
1849  do ij = 1, icells
1850  i = indxi(ij)
1851  j = indxj(ij)
1852 
1853  ! update fluxes that depend on Tsf
1854  flwoutn(i,j) = flwoutn(i,j) + dtsf(i,j) * dflwout_dt(i,j)
1855  fsensn(i,j) = fsensn(i,j) + dtsf(i,j) * dfsens_dt(i,j)
1856  flatn(i,j) = flatn(i,j) + dtsf(i,j) * dflat_dt(i,j)
1857 
1858  ! absorbed shortwave flux for coupler
1859  fswabsn(i,j) = fswsfc(i,j) + fswint(i,j) + fswthrun(i,j)
1860 
1861  enddo ! ij
1862 
1863  end subroutine temperature_changes
1864 
1865 !=======================================================================
1866 !BOP
1867 !
1868 ! !ROUTINE: conductivity - compute ice thermal conductivity
1869 !
1870 ! !DESCRIPTION:
1871 !
1872 ! Compute thermal conductivity at interfaces (held fixed during
1873 ! the subsequent iteration).
1874 !
1875 ! NOTE: Ice conductivity must be >= kimin
1876 !
1877 ! !REVISION HISTORY:
1878 !
1879 ! authors William H. Lipscomb, LANL
1880 ! C. M. Bitz, UW
1881 !
1882 ! !INTERFACE:
1883 !
1884  subroutine conductivity &
1885  (icells, indxi, indxj, &
1886  hlyr, hsn, tin, tbot, khi, khs)
1887 !
1888 ! !USES:
1889 !
1890 ! !INPUT/OUTPUT PARAMETERS:
1891 !
1892  integer (kind=int_kind), intent(in) :: &
1893  icells ! number of cells with aicen > puny
1894 
1895  integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), &
1896  intent(in) :: &
1897  indxi, indxj ! compressed indices for cells with aicen > puny
1898 
1899  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(in) :: &
1900  hlyr & ! ice layer thickness
1901  , hsn & ! snow layer thickness
1902  , tbot ! ice bottom surface temperature (deg C)
1903 
1904  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nilyr), &
1905  intent(in) :: &
1906  tin ! internal ice layer temperatures
1907 
1908  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nilyr+1), &
1909  intent(out) :: &
1910  khi ! ki / hlyr
1911 
1912  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(out) :: &
1913  khs ! ksno / hsn
1914 !
1915 !EOP
1916 !
1917  real (kind=dbl_kind), parameter :: &
1918  betak = 0.13_dbl_kind &! constant in formula for k (W m-1 ppt-1)
1919  , kimin = 0.10_dbl_kind ! min conductivity of saline ice (W m-1 deg-1)
1920 
1921  integer (kind=int_kind) :: &
1922  i, j & ! horizontal indices
1923  , ij & ! horizontal index, combines i and j loops
1924  , k ! ice layer index
1925 
1926  real (kind=dbl_kind) :: &
1927  ki ! thermal cond of sea ice (W m-1 deg-1)
1928 
1929  khi(:,:,:) = c0i
1930  khs(:,:) = c0i
1931 
1932 !DIR$ CONCURRENT !Cray
1933 !cdir nodep !NEC
1934 !ocl novrec !Fujitsu
1935  do ij = 1, icells
1936  i = indxi(ij)
1937  j = indxj(ij)
1938 
1939  ! top surface of top ice layer
1940  ki = kice + betak*salin(1) / min(-puny,tin(i,j,1))
1941  ki = max(ki, kimin)
1942  khi(i,j,1) = ki / hlyr(i,j)
1943 
1944  ! bottom surface of bottom ice layer
1945  ki = kice + betak*salin(nilyr+1) / min(-puny,tbot(i,j))
1946  ki = max(ki, kimin)
1947  khi(i,j,nilyr+1) = ki / hlyr(i,j)
1948 
1949  ! if snow is present: top snow surface, snow-ice interface
1950  if (hsn(i,j) > hsnomin) then
1951  khs(i,j) = ksno / hsn(i,j)
1952  khi(i,j,1) = c2i*khi(i,j,1)*khs(i,j) / &
1953  (khi(i,j,1) + khs(i,j))
1954  endif
1955 
1956  enddo ! ij
1957 
1958  ! interior ice interfaces
1959  do k = 2, nilyr
1960 !DIR$ CONCURRENT !Cray
1961 !cdir nodep !NEC
1962 !ocl novrec !Fujitsu
1963  do ij = 1, icells
1964  i = indxi(ij)
1965  j = indxj(ij)
1966 
1967  ki = kice + betak*p5*(salin(k-1)+salin(k)) / &
1968  min(-puny, p5*(tin(i,j,k-1)+tin(i,j,k)))
1969  ki = max(ki, kimin)
1970  khi(i,j,k) = ki / hlyr(i,j)
1971 
1972  enddo ! ij
1973  enddo ! k
1974 
1975  end subroutine conductivity
1976 
1977 !=======================================================================
1978 !BOP
1979 !
1980 ! !ROUTINE: absorbed_solar - shortwave radiation absorbed by ice, ocean
1981 !
1982 ! !DESCRIPTION:
1983 !
1984 ! Compute solar radiation absorbed in ice and penetrating to ocean
1985 !
1986 ! !REVISION HISTORY:
1987 !
1988 ! authors William H. Lipscomb, LANL
1989 ! C. M. Bitz, UW
1990 !
1991 ! !INTERFACE:
1992 !
1993  subroutine absorbed_solar &
1994  (ni, icells, indxi, indxj, &
1995  hlyr, hsn, fswsfc, fswint, fswthrun, iabs)
1996 !
1997 ! !USES:
1998 !
1999  use ice_albedo
2000 !
2001 ! !INPUT/OUTPUT PARAMETERS:
2002 !
2003  integer (kind=int_kind), intent(in) :: &
2004  ni & ! thickness category index
2005  , icells ! number of cells with aicen > puny
2006 
2007  integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), &
2008  intent(in) :: &
2009  indxi, indxj ! compressed indices for cells with aicen > puny
2010 
2011  real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(in) :: &
2012  hlyr & ! ice layer thickness
2013  , hsn ! snow thickness (m)
2014 
2015  real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(inout) :: &
2016  fswsfc & ! SW absorbed at ice/snow surface (W m-2)
2017  , fswint & ! SW absorbed in ice interior, below surface (W m-2)
2018  , fswthrun ! SW through ice to ocean (W m-2)
2019 
2020  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nilyr), &
2021  intent(inout) :: &
2022  iabs ! SW absorbed in particular layer (W m-2)
2023 !
2024 !EOP
2025 !
2026  real (kind=dbl_kind), parameter :: &
2027  i0vis = 0.70_dbl_kind ! fraction of penetrating solar rad (visible)
2028 
2029  integer (kind=int_kind) :: &
2030  i, j & ! horizontal indices
2031  , ij & ! horizontal index, combines i and j loops
2032  , k ! ice layer index
2033 
2034  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi) :: &
2035  fswpen & ! SW penetrating beneath surface (W m-2)
2036  , trantop & ! transmitted frac of penetrating SW at layer top
2037  , tranbot ! transmitted frac of penetrating SW at layer bot
2038 
2039  real (kind=dbl_kind) :: &
2040  swabs & ! net SW down at surface (W m-2)
2041  , swabsv & ! swabs in vis (wvlngth < 700nm) (W/m^2)
2042  , swabsi & ! swabs in nir (wvlngth > 700nm) (W/m^2)
2043  , frsnow ! fractional snow coverage
2044 
2045 
2046  fswpen(:,:) = c0i
2047  trantop(:,:) = c0i
2048  tranbot(:,:) = c0i
2049 
2050  do ij = 1, icells
2051  i = indxi(ij)
2052  j = indxj(ij)
2053 
2054  !-----------------------------------------------------------------
2055  ! Fractional snow cover
2056  !-----------------------------------------------------------------
2057 
2058  if (hsn(i,j) > puny) then
2059  frsnow = hsn(i,j) / (hsn(i,j) + snowpatch)
2060  else
2061  frsnow = c0i
2062  endif
2063 
2064  !-----------------------------------------------------------------
2065  ! Shortwave flux absorbed at surface, absorbed internally,
2066  ! and penetrating to mixed layer.
2067  ! This parameterization assumes that all IR is absorbed at the
2068  ! surface; only visible is absorbed in the ice interior or
2069  ! transmitted to the ocean.
2070  !-----------------------------------------------------------------
2071 
2072  swabsv = swvdr(i,j)*(c1i-alvdrn(i,j,ni)) &
2073  + swvdf(i,j)*(c1i-alvdfn(i,j,ni))
2074  swabsi = swidr(i,j)*(c1i-alidrn(i,j,ni)) &
2075  + swidf(i,j)*(c1i-alidfn(i,j,ni))
2076  swabs = swabsv + swabsi
2077 
2078  fswpen(i,j) = swabsv * (c1i-frsnow) * i0vis
2079 !c & + swabsi * (c1i-frsnow) * i0nir ! i0nir = 0
2080  fswsfc(i,j) = swabs - fswpen(i,j)
2081 
2082  trantop(i,j) = c1i ! transmittance at top of ice
2083 
2084  enddo ! ij
2085 
2086  !-----------------------------------------------------------------
2087  ! penetrating SW absorbed in each ice layer
2088  !-----------------------------------------------------------------
2089 
2090  do k = 1, nilyr
2091 !DIR$ CONCURRENT !Cray
2092 !cdir nodep !NEC
2093 !ocl novrec !Fujitsu
2094  do ij = 1, icells
2095  i = indxi(ij)
2096  j = indxj(ij)
2097 
2098  tranbot(i,j) = exp(-kappav*hlyr(i,j)*real(k,kind=dbl_kind))
2099  iabs(i,j,k) = fswpen(i,j) * (trantop(i,j)-tranbot(i,j))
2100 
2101  ! bottom of layer k = top of layer k+1
2102  trantop(i,j) = tranbot(i,j)
2103 
2104  enddo ! ij
2105  enddo ! k
2106 
2107 !DIR$ CONCURRENT !Cray
2108 !cdir nodep !NEC
2109 !ocl novrec !Fujitsu
2110  do ij = 1, icells
2111  i = indxi(ij)
2112  j = indxj(ij)
2113 
2114  ! SW penetrating thru ice into ocean
2115  fswthrun(i,j) = fswpen(i,j) * tranbot(i,j)
2116 
2117  ! SW absorbed in ice interior
2118  fswint(i,j) = fswpen(i,j) - fswthrun(i,j)
2119 
2120  enddo ! ij
2121 
2122  end subroutine absorbed_solar
2123 
2124 !=======================================================================
2125 !BOP
2126 !
2127 ! !ROUTINE: surface_fluxes - surface radiative and turbulent fluxes
2128 !
2129 ! !DESCRIPTION:
2130 !
2131 ! Compute radiative and turbulent fluxes and their derivatives
2132 ! with respect to Tsf.
2133 !
2134 ! !REVISION HISTORY:
2135 !
2136 ! authors William H. Lipscomb, LANL
2137 ! C. M. Bitz, UW
2138 !
2139 ! !INTERFACE:
2140 !
2141  subroutine surface_fluxes &
2142  (isolve, indxii, indxjj, &
2143  tsf, fswsfc, &
2144  flwoutn, fsensn, flatn, fsurf, &
2145  dflwout_dt, dfsens_dt, dflat_dt, dfsurf_dt)
2146 !
2147 ! !USES:
2148 !
2149 ! !INPUT/OUTPUT PARAMETERS:
2150 !
2151  integer (kind=int_kind), intent(in) :: &
2152  isolve ! number of cells with temps not converged
2153 
2154  integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), &
2155  intent(in) :: &
2156  indxii, indxjj ! compressed indices for cells not converged
2157 
2158  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(in) :: &
2159  tsf & ! ice/snow surface temperature, Tsfcn
2160  , fswsfc ! SW absorbed at ice/snow surface (W m-2)
2161 
2162  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout):: &
2163  fsensn & ! surface downward sensible heat (W m-2)
2164  , flatn & ! surface downward latent heat (W m-2)
2165  , flwoutn & ! upward LW at surface (W m-2)
2166  , fsurf & ! net flux to top surface, not including fcondtop
2167  , dfsens_dt & ! deriv of fsens wrt Tsf (W m-2 deg-1)
2168  , dflat_dt & ! deriv of flat wrt Tsf (W m-2 deg-1)
2169  , dflwout_dt & ! deriv of flwout wrt Tsf (W m-2 deg-1)
2170  , dfsurf_dt ! derivative of fsurf wrt Tsf
2171 !
2172 !EOP
2173 !
2174  integer (kind=int_kind) :: &
2175  i, j & ! horizontal indices
2176  , ij & ! horizontal index, combines i and j loops
2177  , k ! ice layer index
2178 
2179  real (kind=dbl_kind) :: &
2180  tsfk & ! ice/snow surface temperature (K)
2181  , qsfc & ! saturated surface specific humidity (kg/kg)
2182  , dqsfcdt & ! derivative of Qsfc wrt surface temperature
2183  , qsat & ! the saturation humidity of air (kg/m^3)
2184  , flwdabs & ! downward longwave absorbed heat flx (W/m^2)
2185  , tmpvar ! 1/TsfK
2186 
2187 !DIR$ CONCURRENT !Cray
2188 !cdir nodep !NEC
2189 !ocl novrec !Fujitsu
2190  do ij = 1, isolve
2191  i = indxii(ij) ! NOTE: not indxi and indxj
2192  j = indxjj(ij)
2193 
2194 
2195  ! ice surface temperature in Kelvin
2196  tsfk = tsf(i,j) + tffresh
2197  tmpvar = c1i/tsfk
2198 
2199  ! saturation humidity
2200  qsat = qqqice * exp(-tttice*tmpvar)
2201  qsfc = qsat / rhoa(i,j)
2202  dqsfcdt = tttice * tmpvar*tmpvar * qsfc
2203 
2204  ! longwave radiative flux
2205  flwdabs = emissivity * flw(i,j)
2206  flwoutn(i,j) = -emissivity * stefan_boltzmann * tsfk**4
2207 
2208  ! downward latent and sensible heat fluxes
2209  fsensn(i,j) = shcoef(i,j) * (pott(i,j) - tsfk)
2210  flatn(i,j) = lhcoef(i,j) * (qa(i,j) - qsfc)
2211 
2212  ! derivatives wrt surface temp
2213  dflwout_dt(i,j) = - emissivity*stefan_boltzmann * c4i*tsfk**3
2214  dfsens_dt(i,j) = - shcoef(i,j)
2215  dflat_dt(i,j) = - lhcoef(i,j) * dqsfcdt
2216 
2217  fsurf(i,j) = fswsfc(i,j) + flwdabs + flwoutn(i,j) &
2218  + fsensn(i,j) + flatn(i,j)
2219  dfsurf_dt(i,j) = dflwout_dt(i,j) &
2220  + dfsens_dt(i,j) + dflat_dt(i,j)
2221  enddo
2222 
2223  end subroutine surface_fluxes
2224 
2225 !=======================================================================
2226 !BOP
2227 !
2228 ! !ROUTINE: tridiag_solver - tridiagonal matrix solver
2229 !
2230 ! !DESCRIPTION:
2231 !
2232 ! Tridiagonal matrix solver--used to solve the implicit vertical heat
2233 ! equation in ice and snow
2234 !
2235 ! !REVISION HISTORY:
2236 !
2237 ! authors William H. Lipscomb, LANL
2238 ! C. M. Bitz, UW
2239 !
2240 ! !INTERFACE:
2241 !
2242  subroutine tridiag_solver &
2243  (isolve, indxii, indxjj, &
2244  nmat, sbdiag, diag, spdiag, rhs, xout)
2245 !
2246 ! !USES:
2247 !
2248 ! !INPUT/OUTPUT PARAMETERS:
2249 !
2250  integer (kind=int_kind), intent(in) :: &
2251  isolve ! number of cells with temps not converged
2252 
2253  integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)),&
2254  intent(in) :: &
2255  indxii, indxjj ! compressed indices for cells not converged
2256 
2257  integer (kind=int_kind), intent(in) :: &
2258  nmat ! matrix dimension
2259 
2260  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nmat), &
2261  intent(in) :: &
2262  diag & ! diagonal matrix elements
2263  , sbdiag & ! sub-diagonal matrix elements
2264  , spdiag & ! super-diagonal matrix elements
2265  , rhs ! rhs of tri-diagonal matrix eqn.
2266 
2267  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nmat), &
2268  intent(out) :: &
2269  xout ! solution vector
2270 !
2271 !EOP
2272 !
2273  integer (kind=int_kind) :: &
2274  i, j & ! horizontal indices
2275  , ij & ! horizontal index, combines i and j loops
2276  , k ! row counter
2277 
2278  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi) :: &
2279  wbeta ! temporary matrix variable
2280 
2281  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nmat) :: &
2282  wgamma ! temporary matrix variable
2283 
2284 
2285 !DIR$ CONCURRENT !Cray
2286 !cdir nodep !NEC
2287 !ocl novrec !Fujitsu
2288  do ij = 1, isolve
2289  i = indxii(ij)
2290  j = indxjj(ij)
2291 
2292  wbeta(i,j) = diag(i,j,1)
2293  xout(i,j,1) = rhs(i,j,1) / wbeta(i,j)
2294 
2295  enddo ! ij
2296 
2297  do k = 2, nmat
2298 !DIR$ CONCURRENT !Cray
2299 !cdir nodep !NEC
2300 !ocl novrec !Fujitsu
2301  do ij = 1, isolve
2302  i = indxii(ij)
2303  j = indxjj(ij)
2304 
2305  wgamma(i,j,k) = spdiag(i,j,k-1) / wbeta(i,j)
2306  wbeta(i,j) = diag(i,j,k) - sbdiag(i,j,k)*wgamma(i,j,k)
2307  xout(i,j,k) = (rhs(i,j,k) - sbdiag(i,j,k)*xout(i,j,k-1)) &
2308  / wbeta(i,j)
2309 
2310  enddo ! ij
2311  enddo ! k
2312 
2313  do k = nmat-1, 1, -1
2314 !DIR$ CONCURRENT !Cray
2315 !cdir nodep !NEC
2316 !ocl novrec !Fujitsu
2317  do ij = 1, isolve
2318  i = indxii(ij)
2319  j = indxjj(ij)
2320 
2321  xout(i,j,k) = xout(i,j,k) - wgamma(i,j,k+1)*xout(i,j,k+1)
2322 
2323  enddo ! ij
2324  enddo ! k
2325 
2326  end subroutine tridiag_solver
2327 
2328 !=======================================================================
2329 !BOP
2330 !
2331 ! !ROUTINE: thickness changes - top and bottom growth/melting
2332 !
2333 ! !DESCRIPTION:
2334 !
2335 ! Compute growth and/or melting at the top and bottom surfaces.
2336 !
2337 ! !REVISION HISTORY:
2338 !
2339 ! authors William H. Lipscomb, LANL
2340 ! C. M. Bitz, UW
2341 !
2342 ! !INTERFACE:
2343 !
2344  subroutine thickness_changes &
2345  (ni, icells, indxi, indxj, &
2346  hin, hlyr, hsn, qin, qsn, &
2347  mvap, tbot, fbot, flatn, fhnetn, &
2348  fsurf, fcondtop, fcondbot, efinal)
2349 !
2350 ! !USES:
2351 !
2352 ! !INPUT/OUTPUT PARAMETERS:
2353 !
2354  integer (kind=int_kind), intent(in) :: &
2355  ni & ! thickness category index
2356  , icells ! number of cells with aicen > puny
2357 
2358  integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), &
2359  intent(in) :: &
2360  indxi, indxj ! compressed indices for cells with aicen > puny
2361 
2362  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(in) :: &
2363  fbot & ! ice-ocean heat flux at bottom surface (W/m^2)
2364  , tbot & ! ice bottom surface temperature (deg C)
2365  , flatn & ! surface downward latent heat (W m-2)
2366  , fsurf & ! net flux to top surface, not including fcondtop
2367  , fcondtop & ! downward cond flux at top surface (W m-2)
2368  , fcondbot ! downward cond flux at bottom surface (W m-2)
2369 
2370  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nilyr), &
2371  intent(inout) :: &
2372  qin ! ice layer enthalpy (J m-3)
2373 
2374  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout):: &
2375  fhnetn & ! fbot, corrected for any surplus energy
2376  , hlyr & ! ice layer thickness
2377  , hsn & ! snow thickness (m)
2378  , qsn & ! snow enthalpy
2379  , efinal & ! final energy of melting (J m-2)
2380  , mvap ! ice/snow mass sublimated/condensed (kg m-2)
2381 
2382  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(out) :: &
2383  hin ! total ice thickness
2384 !
2385 !EOP
2386 !
2387  real (kind=dbl_kind), parameter :: &
2388  qbotmax = -p5*rhoi*lfresh ! max enthalpy of ice growing at bottom
2389 
2390  integer (kind=int_kind) :: &
2391  i, j & ! horizontal indices
2392  , ij & ! horizontal index, combines i and j loops
2393  , k, kold ! ice layer indices
2394 
2395  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi) :: &
2396  esub &! energy for sublimation, > 0 (J m-2)
2397  , etop_mlt &! energy for top melting, > 0 (J m-2)
2398  , ebot_mlt &! energy for bottom melting, > 0 (J m-2)
2399  , rhlyr ! reciprocal ice layer thickness
2400 
2401  real (kind=dbl_kind) :: &
2402  dhi & ! change in ice thickness
2403  , dhs & ! change in snow thickness
2404  , ti & ! ice temperature
2405  , ts & ! snow temperature
2406  , econ & ! energy for condensation, < 0 (J m-2)
2407  , ebot_gro & ! energy for bottom growth, < 0 (J m-2)
2408  , qbot & ! enthalpy of ice growing at bottom surface (J m-3)
2409  , qsub & ! energy/unit volume to sublimate ice/snow (J m-3)
2410  , hqtot & ! sum of h*q for two layers
2411  , w1 & ! temporary variable
2412  , hovlp ! overlap between old and new layers (m)
2413 
2414  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nilyr) :: &
2415  hnew & ! new layer thickness
2416  , hq ! h * q for a layer
2417 
2418  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,0:nilyr) :: &
2419  zold & ! depth of layer boundaries (m)
2420  , znew ! adjusted depths, with equal hlyr (m)
2421 
2422  !-----------------------------------------------------------------
2423  ! Initialize ice layer thicknesses
2424  !-----------------------------------------------------------------
2425 
2426  do k = 1, nilyr
2427  do ij = 1, icells
2428  i = indxi(ij)
2429  j = indxj(ij)
2430  hnew(i,j,k) = hlyr(i,j)
2431  enddo ! if
2432  enddo
2433 
2434  !-----------------------------------------------------------------
2435  ! For l_brine = false (fresh ice), check for temperatures > 0.
2436  ! Melt ice or snow as needed to bring temperatures back to 0.
2437  ! For l_brine = true, this should not be necessary.
2438  !-----------------------------------------------------------------
2439 
2440  if (.not. l_brine) then
2441 !DIR$ CONCURRENT !Cray
2442 !cdir nodep !NEC
2443 !ocl novrec !Fujitsu
2444  do ij = 1, icells
2445  i = indxi(ij)
2446  j = indxj(ij)
2447 
2448  ts = (lfresh + qsn(i,j)/rhos) / cp_ice
2449  if (ts > c0i) then
2450  dhs = cp_ice*ts*hsn(i,j) / lfresh
2451  hsn(i,j) = hsn(i,j) - dhs
2452  qsn(i,j) = -rhos*lfresh
2453  endif
2454  enddo
2455 
2456  do k = 1, nilyr
2457 !DIR$ CONCURRENT !Cray
2458 !cdir nodep !NEC
2459 !ocl novrec !Fujitsu
2460  do ij = 1, icells
2461  i = indxi(ij)
2462  j = indxj(ij)
2463 
2464  ti = (lfresh + qin(i,j,k)/rhoi) / cp_ice
2465  if (ti > c0i) then
2466  dhi = cp_ice*ti*hnew(i,j,k) / lfresh
2467  hnew(i,j,k) = hnew(i,j,k) - dhi
2468  qin(i,j,k) = -rhoi*lfresh
2469  endif
2470  enddo ! ij
2471  enddo ! k
2472 
2473  endif ! .not. l_brine
2474 
2475  !-----------------------------------------------------------------
2476  ! Sublimation/condensation of ice/snow
2477  !-----------------------------------------------------------------
2478 
2479 !DIR$ CONCURRENT !Cray
2480 !cdir nodep !NEC
2481 !ocl novrec !Fujitsu
2482  do ij = 1, icells
2483  i = indxi(ij)
2484  j = indxj(ij)
2485 
2486 ! w1 = -flatn(i,j) * dt
2487  w1 = -flatn(i,j) * dtice
2488  esub(i,j) = max(w1, c0i) ! energy for sublimation, > 0
2489  econ = min(w1, c0i) ! energy for condensation, < 0
2490 
2491  !--------------------------------------------------------------
2492  ! Condensation (mvap > 0)
2493  !--------------------------------------------------------------
2494 
2495  if (hsn(i,j) > puny) then ! add snow with enthalpy qsn
2496  dhs = econ / (qsn(i,j) - rhos*lvap) ! econ < 0, dhs > 0
2497  hsn(i,j) = hsn(i,j) + dhs
2498  mvap(i,j) = mvap(i,j) + dhs*rhos
2499  else ! add ice with enthalpy qin(i,j,1)
2500  dhi = econ / (qin(i,j,1) - rhoi*lvap) ! econ < 0, dhi > 0
2501  hnew(i,j,1) = hnew(i,j,1) + dhi
2502  mvap(i,j) = mvap(i,j) + dhi*rhoi
2503  endif
2504 
2505  !--------------------------------------------------------------
2506  ! Sublimation of snow (mvap < 0)
2507  !--------------------------------------------------------------
2508  qsub = qsn(i,j) - rhos*lvap ! qsub < 0
2509  dhs = max(-hsn(i,j), esub(i,j)/qsub) ! esub > 0, dhs < 0
2510  hsn(i,j) = hsn(i,j) + dhs
2511  esub(i,j) = esub(i,j) - dhs*qsub
2512  esub(i,j) = max(esub(i,j), c0i) ! in case of roundoff error
2513  mvap(i,j) = mvap(i,j) + dhs*rhos
2514  enddo ! ij
2515 
2516  !--------------------------------------------------------------
2517  ! Sublimation of ice (mvap < 0)
2518  !--------------------------------------------------------------
2519 
2520  do k = 1, nilyr
2521 !DIR$ CONCURRENT !Cray
2522 !cdir nodep !NEC
2523 !ocl novrec !Fujitsu
2524  do ij = 1, icells
2525  i = indxi(ij)
2526  j = indxj(ij)
2527  qsub = qin(i,j,k) - rhoi*lvap ! qsub < 0
2528  dhi = max(-hnew(i,j,k), esub(i,j)/qsub) ! esub < 0, dhi < 0
2529  hnew(i,j,k) = hnew(i,j,k) + dhi
2530  esub(i,j) = esub(i,j) - dhi*qsub
2531  esub(i,j) = max(esub(i,j), c0i)
2532  mvap(i,j) = mvap(i,j) + dhi*rhoi
2533  enddo ! ij
2534  enddo ! k
2535 
2536  !-----------------------------------------------------------------
2537  ! Top melt
2538  ! (There is no top growth because there is no liquid water to
2539  ! freeze at the top.)
2540  !-----------------------------------------------------------------
2541 
2542  !--------------------------------------------------------------
2543  ! Melt snow
2544  !--------------------------------------------------------------
2545 !DIR$ CONCURRENT !Cray
2546 !cdir nodep !NEC
2547 !ocl novrec !Fujitsu
2548  do ij = 1, icells
2549  i = indxi(ij)
2550  j = indxj(ij)
2551 
2552 ! w1 = (fsurf(i,j) - fcondtop(i,j)) * dt
2553  w1 = (fsurf(i,j) - fcondtop(i,j)) * dtice
2554  etop_mlt(i,j) = max(w1, c0i) ! etop_mlt > 0
2555  dhs = max(-hsn(i,j), etop_mlt(i,j)/qsn(i,j)) ! qsn < 0, dhs < 0
2556  hsn(i,j) = hsn(i,j) + dhs
2557  etop_mlt(i,j) = etop_mlt(i,j) - dhs*qsn(i,j)
2558  etop_mlt(i,j) = max(etop_mlt(i,j), c0i) ! in case of roundoff error
2559 
2560  ! history diagnostics
2561  if (dhs < -puny .and. mlt_onset(i,j) < puny) &
2562  mlt_onset(i,j) = yday
2563 
2564  enddo ! ij
2565 
2566  !--------------------------------------------------------------
2567  ! Melt ice
2568  !--------------------------------------------------------------
2569 
2570  do k = 1, nilyr
2571 !DIR$ CONCURRENT !Cray
2572 !cdir nodep !NEC
2573 !ocl novrec !Fujitsu
2574  do ij = 1, icells
2575  i = indxi(ij)
2576  j = indxj(ij)
2577 
2578  dhi = max(-hnew(i,j,k), etop_mlt(i,j)/qin(i,j,k))
2579  hnew(i,j,k) = hnew(i,j,k) + dhi ! qin < 0, dhi < 0
2580  etop_mlt(i,j) = etop_mlt(i,j) - dhi*qin(i,j,k)
2581  etop_mlt(i,j) = max(etop_mlt(i,j), c0i)
2582 
2583  ! history diagnostics
2584  if (dhi < -puny .and. mlt_onset(i,j) < puny) &
2585  mlt_onset(i,j) = yday
2586  meltt(i,j) = meltt(i,j) - dhi*aicen(i,j,ni)
2587 
2588  enddo ! ij
2589  enddo ! k
2590 
2591  !----------------------------------------------------------------
2592  ! Bottom growth and melt
2593  !----------------------------------------------------------------
2594 
2595 !DIR$ CONCURRENT !Cray
2596 !cdir nodep !NEC
2597 !ocl novrec !Fujitsu
2598  do ij = 1, icells
2599  i = indxi(ij)
2600  j = indxj(ij)
2601 
2602 ! w1 = (fcondbot(i,j) - fbot(i,j)) * dt
2603  w1 = (fcondbot(i,j) - fbot(i,j)) * dtice
2604  ebot_mlt(i,j) = max(w1, c0i) ! ebot_mlt > 0
2605  ebot_gro = min(w1, c0i) ! ebot_gro < 0
2606 
2607  !--------------------------------------------------------------
2608  ! Grow ice
2609  !--------------------------------------------------------------
2610 
2611  ! enthalpy of new ice growing at bottom surface
2612  qbot = -rhoi * (cp_ice * (tmlt(nilyr+1)-tbot(i,j)) &
2613  + lfresh * (c1i-tmlt(nilyr+1)/tbot(i,j)) &
2614  - cp_ocn * tmlt(nilyr+1))
2615  qbot = min(qbot, qbotmax) ! in case Tbot is close to Tmlt
2616  dhi = ebot_gro / qbot ! dhi > 0
2617 
2618  hqtot = hnew(i,j,nilyr)*qin(i,j,nilyr) + dhi*qbot
2619  hnew(i,j,nilyr) = hnew(i,j,nilyr) + dhi
2620 
2621  if (hnew(i,j,nilyr) > puny) &
2622  qin(i,j,nilyr) = hqtot / hnew(i,j,nilyr)
2623 
2624  ! history diagnostics
2625  congel(i,j) = congel(i,j) + dhi*aicen(i,j,ni)
2626  if (dhi > puny .and. frz_onset(i,j) < puny) &
2627  frz_onset(i,j) = yday
2628 
2629  enddo ! ij
2630 
2631 
2632  !--------------------------------------------------------------
2633  ! Melt ice
2634  !--------------------------------------------------------------
2635 
2636  do k = nilyr, 1, -1
2637 !DIR$ CONCURRENT !Cray
2638 !cdir nodep !NEC
2639 !ocl novrec !Fujitsu
2640  do ij = 1, icells
2641  i = indxi(ij)
2642  j = indxj(ij)
2643 
2644  dhi = max(-hnew(i,j,k), ebot_mlt(i,j)/qin(i,j,k))
2645  hnew(i,j,k) = hnew(i,j,k) + dhi ! qin < 0, dhi < 0
2646  ebot_mlt(i,j) = ebot_mlt(i,j) - dhi*qin(i,j,k)
2647  ebot_mlt(i,j) = max(ebot_mlt(i,j), c0i)
2648 
2649  ! history diagnostics
2650  meltb(i,j) = meltb(i,j) - dhi*aicen(i,j,ni)
2651 
2652  enddo ! ij
2653  enddo ! k
2654 
2655  !--------------------------------------------------------------
2656  ! Melt snow (only if all the ice has melted)
2657  !--------------------------------------------------------------
2658 
2659 !DIR$ CONCURRENT !Cray
2660 !cdir nodep !NEC
2661 !ocl novrec !Fujitsu
2662  do ij = 1, icells
2663  i = indxi(ij)
2664  j = indxj(ij)
2665 
2666  dhs = max(-hsn(i,j), ebot_mlt(i,j)/qsn(i,j))
2667  hsn(i,j) = hsn(i,j) + dhs ! qsn < 0, dhs < 0
2668  ebot_mlt(i,j) = ebot_mlt(i,j) - dhs*qsn(i,j)
2669  ebot_mlt(i,j) = max(ebot_mlt(i,j), c0i)
2670 
2671  ! snow-to-ice conversion (very rare)
2672  if (hsn(i,j) > puny .and. ebot_mlt(i,j) > puny*lfresh) then
2673  hnew(i,j,1) = hsn(i,j) * rhos/rhoi ! reduce thickness
2674  qin(i,j,1) = qsn(i,j) * rhoi/rhos ! increase q to conserve energy
2675  hsn(i,j) = c0i
2676  endif
2677 
2678  enddo ! ij
2679 
2680  !-----------------------------------------------------------------
2681  ! Give the ocean any energy left over after melting
2682  !-----------------------------------------------------------------
2683 
2684 !DIR$ CONCURRENT !Cray
2685 !cdir nodep !NEC
2686 !ocl novrec !Fujitsu
2687  do ij = 1, icells
2688  i = indxi(ij)
2689  j = indxj(ij)
2690 
2691  fhnetn(i,j) = fhnetn(i,j) &
2692  + (esub(i,j) + etop_mlt(i,j) + ebot_mlt(i,j))/dtice
2693 ! + (esub(i,j) + etop_mlt(i,j) + ebot_mlt(i,j))/dt
2694 
2695  enddo ! ij
2696 
2697 !---!-------------------------------------------------------------------
2698 !---! Repartition the ice into equal-thickness layers, conserving energy.
2699 !---!-------------------------------------------------------------------
2700 
2701  !-----------------------------------------------------------------
2702  ! Initialize the new ice thickness.
2703  ! Initialize the final energy: snow energy plus energy of
2704  ! sublimated/condensed ice.
2705  !-----------------------------------------------------------------
2706 
2707  do ij = 1, icells
2708  i = indxi(ij)
2709  j = indxj(ij)
2710  hin(i,j) = c0i
2711  efinal(i,j) = hsn(i,j)*qsn(i,j) - mvap(i,j)*lvap
2712  enddo
2713 
2714  !-----------------------------------------------------------------
2715  ! Compute the new ice thickness.
2716  ! Initialize h*q for new layers.
2717  !-----------------------------------------------------------------
2718 
2719  do k = 1, nilyr
2720 !DIR$ CONCURRENT !Cray
2721 !cdir nodep !NEC
2722 !ocl novrec !Fujitsu
2723  do ij = 1, icells
2724  i = indxi(ij)
2725  j = indxj(ij)
2726  hin(i,j) = hin(i,j) + hnew(i,j,k)
2727  hq(i,j,k) = c0i
2728  enddo ! ij
2729  enddo ! k
2730 
2731  !-----------------------------------------------------------------
2732  ! Make sure hin > puny.
2733  ! Compute depths zold of old layers (unequal thicknesses).
2734  ! Compute depths znew of new layers (all the same thickness).
2735  !-----------------------------------------------------------------
2736 
2737 !DIR$ CONCURRENT !Cray
2738 !cdir nodep !NEC
2739 !ocl novrec !Fujitsu
2740  do ij = 1, icells
2741  i = indxi(ij)
2742  j = indxj(ij)
2743 
2744  if (hin(i,j) > puny) then
2745  hlyr(i,j) = hin(i,j) / real(nilyr,kind=dbl_kind)
2746  rhlyr(i,j) = c1i / hlyr(i,j)
2747  else
2748  hin(i,j) = c0i
2749  hlyr(i,j) = c0i
2750  rhlyr(i,j) = c0i
2751  endif
2752 
2753  zold(i,j,0) = c0i
2754  znew(i,j,0) = c0i
2755 
2756  zold(i,j,nilyr) = hin(i,j)
2757  znew(i,j,nilyr) = hin(i,j)
2758 
2759  enddo ! ij
2760 
2761  do k = 1, nilyr-1
2762  do ij = 1, icells
2763  i = indxi(ij)
2764  j = indxj(ij)
2765 
2766  zold(i,j,k) = zold(i,j,k-1) + hnew(i,j,k) ! old unequal layers
2767  znew(i,j,k) = znew(i,j,k-1) + hlyr(i,j) ! new equal layers
2768 
2769  end do ! ij
2770  enddo ! k
2771 
2772  !-----------------------------------------------------------------
2773  ! Compute h*q for new layer (k) given overlap with old layers (kold)
2774  !-----------------------------------------------------------------
2775 
2776  do k = 1, nilyr
2777  do kold = 1, nilyr
2778 !DIR$ CONCURRENT !Cray
2779 !cdir nodep !NEC
2780 !ocl novrec !Fujitsu
2781  do ij = 1, icells
2782  i = indxi(ij)
2783  j = indxj(ij)
2784 
2785  hovlp = min(zold(i,j,kold), znew(i,j,k)) &
2786  - max(zold(i,j,kold-1), znew(i,j,k-1))
2787  hovlp = max(hovlp, c0i)
2788 
2789  hq(i,j,k) = hq(i,j,k) + hovlp*qin(i,j,kold)
2790 
2791  enddo ! ij
2792  enddo ! kold
2793  enddo ! k
2794 
2795  !-----------------------------------------------------------------
2796  ! Compute new values of qin and increment ice-snow energy.
2797  ! Note: Have to finish previous loop before updating qin
2798  ! in this loop.
2799  !-----------------------------------------------------------------
2800 
2801  do k = 1, nilyr
2802 !DIR$ CONCURRENT !Cray
2803 !cdir nodep !NEC
2804 !ocl novrec !Fujitsu
2805  do ij = 1, icells
2806  i = indxi(ij)
2807  j = indxj(ij)
2808 
2809  qin(i,j,k) = hq(i,j,k) * rhlyr(i,j)
2810  efinal(i,j) = efinal(i,j) + hlyr(i,j)*qin(i,j,k)
2811 
2812  enddo ! ij
2813  enddo !
2814 
2815  end subroutine thickness_changes
2816 
2817 !=======================================================================
2818 !BOP
2819 !
2820 ! !ROUTINE: conservation_check_vthermo - energy conservation check
2821 !
2822 ! !DESCRIPTION:
2823 !
2824 ! Check for energy conservation by comparing the change in energy
2825 ! to the net energy input.
2826 !
2827 ! !REVISION HISTORY:
2828 !
2829 ! authors William H. Lipscomb, LANL
2830 ! C. M. Bitz, UW
2831 !
2832 ! !INTERFACE:
2833 !
2834  subroutine conservation_check_vthermo &
2835  (ni, icells, indxi, indxj, &
2836  fsurf, flatn, fhnetn, fswint, einit, efinal)
2837 !
2838 ! !USES:
2839 !
2840 ! use ice_exit
2841 !
2842 ! !INPUT/OUTPUT PARAMETERS:
2843 !
2844  integer (kind=int_kind), intent(in) :: &
2845  ni & ! thickness category index (diagnostic only)
2846  , icells ! number of cells with aicen > puny
2847 
2848  integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), &
2849  intent(in) :: &
2850  indxi, indxj ! compressed indices for cells with aicen > puny
2851 
2852  real (kind=dbl_kind), dimension (ilo:ihi, jlo:jhi), intent(in) :: &
2853  fsurf & ! net flux to top surface, not including fcondtop
2854  , flatn & ! surface downward latent heat (W m-2)
2855  , fhnetn & ! fbot, corrected for any surplus energy
2856  , fswint & ! SW absorbed in ice interior, below surface (W m-2)
2857  , einit & ! initial energy of melting (J m-2)
2858  , efinal ! final energy of melting (J m-2)
2859 !
2860 !EOP
2861 !
2862  integer (kind=int_kind) :: &
2863  i, j & ! horizontal indices
2864  , ij ! horizontal index, combines i and j loops
2865 
2866  real (kind=dbl_kind) :: &
2867  einp & ! energy input during timestep (J m-2)
2868  , ferr ! energy conservation error (W m-2)
2869 
2870  logical (kind=log_kind) :: & ! for vector-friendly error checks
2871  ferr_flag ! flag for energy error, ferr > ferrmax
2872 
2873  ferr_flag = .false. ! ferr <= ferrmax, initialization
2874 
2875 !DIR$ CONCURRENT !Cray
2876 !cdir nodep !NEC
2877 !ocl novrec !Fujitsu
2878  do ij = 1, icells
2879  i = indxi(ij)
2880  j = indxj(ij)
2881 
2882  einp = (fsurf(i,j) - flatn(i,j) + fswint(i,j) - fhnetn(i,j)) &
2883  * dtice
2884 ! * dt
2885 ! ferr = abs(efinal(i,j)-einit(i,j)-einp) / dt
2886  ferr = abs(efinal(i,j)-einit(i,j)-einp) / dtice
2887  if (ferr > ferrmax) then
2888  ferr_flag = .true.
2889  endif
2890  enddo
2891  !----------------------------------------------------------------
2892  ! If energy is not conserved, print diagnostics and exit.
2893  !----------------------------------------------------------------
2894  if (ferr_flag .and. dbg_set(dbg_sbrio)) then
2895  do ij = 1, icells
2896  i = indxi(ij)
2897  j = indxj(ij)
2898 
2899  !-----------------------------------------------------------------
2900  ! Note that fsurf - flat = fsw + flw + fsens; i.e., the latent
2901  ! heat is not included in the energy input, since de is the energy
2902  ! change in the system ice + vapor, and the latent heat lost by
2903  ! the ice is equal to that gained by the vapor.
2904  !-----------------------------------------------------------------
2905 
2906  einp = (fsurf(i,j) - flatn(i,j) + fswint(i,j) - &
2907  fhnetn(i,j)) * dtice
2908 ! fhnetn(i,j)) * dt
2909 ! ferr = abs(efinal(i,j)-einit(i,j)-einp) / dt
2910  ferr = abs(efinal(i,j)-einit(i,j)-einp) / dtice
2911  if (ferr > ferrmax) then
2912 ! write(nu_diag,*) 'Energy error, cat',ni, &
2913 ! 'my_task,i,j',my_task,i,j
2914 ! write(nu_diag,*)'wind(I,j)',wind(i,j),shcoef(i,j), potT(i,j)
2915 ! write(nu_diag,*) 'Flux error (W/m^2) =', ferr
2916 ! write(nu_diag,*) 'Energy error (J) =', ferr*dt
2917 ! write(nu_diag,*) 'Energy error (J) =', ferr*dtice
2918 ! write(nu_diag,*) 'Initial energy =', einit(i,j)
2919 ! write(nu_diag,*) 'Final energy =', efinal(i,j)
2920 ! write(nu_diag,*) 'efinal - einit =', &
2921 ! efinal(i,j)-einit(i,j)
2922 ! write(nu_diag,*) 'Input energy =', einp
2923 ! stoplabel = 'ice state at stop'
2924 ! call print_state(stoplabel,i,j)
2925 ! call abort_ice &
2926 ! ('vertical thermo: energy conservation error')
2927  endif
2928  enddo
2929  endif
2930 
2931  end subroutine conservation_check_vthermo
2932 
2933 !=======================================================================
2934 !BOP
2935 !
2936 ! !ROUTINE: add_new_snow - add new snow
2937 !
2938 ! !DESCRIPTION:
2939 !
2940 ! Add new snow at top surface
2941 !
2942 ! !REVISION HISTORY:
2943 !
2944 ! authors William H. Lipscomb, LANL
2945 ! C. M. Bitz, UW
2946 !
2947 ! !INTERFACE:
2948 !
2949  subroutine add_new_snow &
2950  (ni, icells, indxi, indxj, &
2951  hsn, hsn_new, qsn, tsf)
2952 !
2953 ! !USES:
2954 !
2955 ! !INPUT/OUTPUT PARAMETERS:
2956 !
2957  integer (kind=int_kind), intent(in) :: &
2958  ni & ! thickness category index
2959  , icells ! number of cells with aicen > puny
2960 
2961  integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), &
2962  intent(in) :: &
2963  indxi, indxj ! compressed indices for cells with aicen > puny
2964 
2965  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(in):: &
2966  tsf ! ice/snow surface temperature, Tsfcn
2967 
2968  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(inout):: &
2969  hsn & ! snow thickness (m)
2970  , qsn ! snow enthalpy
2971 
2972  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(out) :: &
2973  hsn_new ! thickness of new snow (m)
2974 !
2975 !EOP
2976 !
2977  real (kind=dbl_kind) :: &
2978  qsnew & ! enthalpy of new snow (J kg-1)
2979  , hstot ! snow thickness including new snow (J kg-1)
2980 
2981  integer (kind=int_kind) :: &
2982  i, j & ! horizontal indices
2983  , ij ! horizontal index, combines i and j loops
2984 
2985 
2986  hsn_new(:,:) = c0i
2987 
2988 !DIR$ CONCURRENT !Cray
2989 !cdir nodep !NEC
2990 !ocl novrec !Fujitsu
2991  do ij = 1, icells
2992  i = indxi(ij)
2993  j = indxj(ij)
2994 
2995  !----------------------------------------------------------------
2996  ! NOTE: If heat flux diagnostics are to work, new snow should
2997  ! have T = 0 (i.e. q = -rhos*Lfresh) and should not be
2998  ! converted to rain.
2999  !----------------------------------------------------------------
3000 
3001  if (fsnow(i,j) > c0i) then
3002 
3003 ! hsn_new(i,j) = fsnow(i,j)/rhos * dt
3004  hsn_new(i,j) = fsnow(i,j)/rhos * dtice
3005  qsnew = -rhos*lfresh
3006  hstot = hsn(i,j) + hsn_new(i,j)
3007  if (hstot > puny) then
3008  qsn(i,j) = (hsn(i,j) *qsn(i,j) &
3009  + hsn_new(i,j)*qsnew ) / hstot
3010  ! avoid roundoff errors if hsn=0
3011  qsn(i,j) = min(qsn(i,j), -rhos*lfresh)
3012  hsn(i,j) = hstot
3013  endif
3014  endif
3015 
3016  enddo ! ij
3017 
3018  end subroutine add_new_snow
3019 
3020 !=======================================================================
3021 !BOP
3022 !
3023 ! !ROUTINE: update_state_vthermo - new state variables
3024 !
3025 ! !DESCRIPTION:
3026 !
3027 ! Given the vertical thermo state variables (hin, hsn, qin,
3028 ! qsn, Tsf), compute the new ice state variables (vicen, vsnon,
3029 ! eicen, esnon, Tsfcn).
3030 ! Zero out state variables if ice has melted entirely.
3031 !
3032 ! !REVISION HISTORY:
3033 !
3034 ! authors William H. Lipscomb, LANL
3035 ! C. M. Bitz, UW
3036 !
3037 ! !INTERFACE:
3038 !
3039  subroutine update_state_vthermo &
3040  (ni, icells, indxi, indxj, &
3041  hin, hsn, qin, qsn, tsf)
3042 !
3043 ! !USES:
3044 !
3045 ! !INPUT/OUTPUT PARAMETERS:
3046 !
3047  integer (kind=int_kind), intent(in) :: &
3048  ni & ! thickness category index
3049  , icells ! number of cells with aicen > puny
3050 
3051  integer (kind=int_kind), dimension(1:(ihi-ilo+1)*(jhi-jlo+1)), &
3052  intent(in) :: &
3053  indxi, indxj ! compressed indices for cells with aicen > puny
3054 
3055 
3056  real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(in) :: &
3057  hin & ! ice thickness (m)
3058  , hsn & ! snow thickness (m)
3059  , qsn & ! snow enthalpy
3060  , tsf ! ice/snow surface temperature, Tsfcn
3061 
3062  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,nilyr), &
3063  intent(in) :: &
3064  qin ! ice layer enthalpy (J m-3)
3065 !
3066 !EOP
3067 !
3068  integer (kind=int_kind) :: &
3069  i, j & ! horizontal indices
3070  , ij & ! horizontal index, combines i and j loops
3071  , k ! ice layer index
3072 
3073 !DIR$ CONCURRENT !Cray
3074 !cdir nodep !NEC
3075 !ocl novrec !Fujitsu
3076  do ij = 1, icells
3077  i = indxi(ij)
3078  j = indxj(ij)
3079 
3080  if (hin(i,j) > c0i) then
3081  ! aicen is already up to date
3082  vicen(i,j,ni) = aicen(i,j,ni) * hin(i,j)
3083  vsnon(i,j,ni) = aicen(i,j,ni) * hsn(i,j)
3084  esnon(i,j,ni) = qsn(i,j) * vsnon(i,j,ni)
3085  tsfcn(i,j,ni) = tsf(i,j)
3086  else ! (hin(i,j) == c0i)
3087  aicen(i,j,ni) = c0i
3088  vicen(i,j,ni) = c0i
3089  vsnon(i,j,ni) = c0i
3090  esnon(i,j,ni) = c0i
3091  tsfcn(i,j,ni) = tf(i,j)
3092  endif
3093 
3094  enddo ! ij
3095 
3096  do k = 1, nilyr
3097 !DIR$ CONCURRENT !Cray
3098 !cdir nodep !NEC
3099 !ocl novrec !Fujitsu
3100  do ij = 1, icells
3101  i = indxi(ij)
3102  j = indxj(ij)
3103 
3104  if (hin(i,j) > c0i) then
3105  eicen(i,j,ilyr1(ni)+k-1) = &
3106  qin(i,j,k) * vicen(i,j,ni)/real(nilyr,kind=dbl_kind)
3107  else
3108  eicen(i,j,ilyr1(ni)+k-1) = c0i
3109  endif
3110 
3111  enddo ! ij
3112  enddo ! k
3113 
3114  end subroutine update_state_vthermo
3115 
3116 !=======================================================================
3117 
3118  end module ice_therm_vertical
3119 
3120 !=======================================================================
ice_flux::tf
real(kind=dbl_kind), dimension(:,:), allocatable, save tf
Definition: ice_flux.f90:91
ice_model_size::nilyr
integer(kind=int_kind), parameter nilyr
Definition: ice_model_size.f90:48
ice_flux::lhcoef
real(kind=dbl_kind), dimension(:,:), allocatable, save lhcoef
Definition: ice_flux.f90:164
ice_constants::depresst
real(kind=dbl_kind), parameter depresst
Definition: ice_constants.f90:118
ice_albedo::alidrn
real(kind=dbl_kind), dimension(:,:,:), allocatable, save alidrn
Definition: ice_albedo.f90:76
mod_utils
Definition: mod_utils.f90:40
ice_state::esnon
real(kind=dbl_kind), dimension(:,:,:), allocatable, target, save esnon
Definition: ice_state.f90:97
ice_constants::tttice
real(kind=dbl_kind), parameter tttice
Definition: ice_constants.f90:134
ice_kinds_mod::dbl_kind
integer, parameter dbl_kind
Definition: ice_kinds_mod.f90:59
mod_utils::dbg_set
logical function dbg_set(vrb)
Definition: mod_utils.f90:182
ice_therm_vertical::hicen_old
real(kind=dbl_kind), dimension(:,:,:), allocatable, save hicen_old
Definition: ice_therm_vertical.f90:84
ice_state
Definition: ice_state.f90:64
ice_therm_vertical::frzmlt_bottom_lateral
subroutine frzmlt_bottom_lateral(Tbot, fbot, rside)
Definition: ice_therm_vertical.f90:477
ice_flux::strocnxt
real(kind=dbl_kind), dimension(:,:), allocatable, save strocnxt
Definition: ice_flux.f90:55
ice_domain::ihi
integer(kind=int_kind) ihi
Definition: ice_domain.f90:101
ice_constants::c0i
real(kind=dbl_kind), parameter c0i
Definition: ice_constants.f90:150
ice_constants::lvap
real(kind=dbl_kind), parameter lvap
Definition: ice_constants.f90:85
ice_constants::c3i
real(kind=dbl_kind), parameter c3i
Definition: ice_constants.f90:150
ice_constants::rhos
real(kind=dbl_kind), parameter rhos
Definition: ice_constants.f90:118
ice_constants::qqqice
real(kind=dbl_kind), parameter qqqice
Definition: ice_constants.f90:134
ice_flux::swvdf
real(kind=dbl_kind), dimension(:,:), allocatable, save swvdf
Definition: ice_flux.f90:91
ice_constants::lfresh
real(kind=dbl_kind), parameter lfresh
Definition: ice_constants.f90:85
ice_constants::ice_ref_salinity
real(kind=dbl_kind), parameter ice_ref_salinity
Definition: ice_constants.f90:85
ice_therm_vertical::init_thermo_vertical
subroutine init_thermo_vertical
Definition: ice_therm_vertical.f90:117
ice_grid
Definition: ice_grid.f90:35
ice_flux::rhoa
real(kind=dbl_kind), dimension(:,:), allocatable, save rhoa
Definition: ice_flux.f90:91
ice_flux::pott
real(kind=dbl_kind), dimension(:,:), allocatable, save pott
Definition: ice_flux.f90:91
ice_therm_vertical::init_vertical_profile
subroutine init_vertical_profile(ni, icells, indxi, indxj, hin, hlyr, hsn, hin_init, hsn_init, qin, Tin, qsn, Tsn, Tsf, einit)
Definition: ice_therm_vertical.f90:763
ice_constants::stefan_boltzmann
real(kind=dbl_kind), parameter stefan_boltzmann
Definition: ice_constants.f90:85
ice_constants::c4i
real(kind=dbl_kind), parameter c4i
Definition: ice_constants.f90:150
ice_albedo::snowpatch
real(kind=dbl_kind), parameter snowpatch
Definition: ice_albedo.f90:59
ice_flux::flw
real(kind=dbl_kind), dimension(:,:), allocatable, save flw
Definition: ice_flux.f90:91
ice_domain::jlo
integer(kind=int_kind) jlo
Definition: ice_domain.f90:101
ice_flux::sst
real(kind=dbl_kind), dimension(:,:), allocatable, save sst
Definition: ice_flux.f90:91
ice_flux::congel
real(kind=dbl_kind), dimension(:,:), allocatable, save congel
Definition: ice_flux.f90:140
ice_flux::init_flux_atm
subroutine init_flux_atm
Definition: ice_flux.f90:182
ice_constants::cp_ice
real(kind=dbl_kind), parameter cp_ice
Definition: ice_constants.f90:85
ice_constants::emissivity
real(kind=dbl_kind), parameter emissivity
Definition: ice_constants.f90:118
ice_therm_vertical::conservation_check_vthermo
subroutine conservation_check_vthermo(ni, icells, indxi, indxj, fsurf, flatn, fhnetn, fswint, einit, efinal)
Definition: ice_therm_vertical.f90:2837
ice_albedo::alidfn
real(kind=dbl_kind), dimension(:,:,:), allocatable, save alidfn
Definition: ice_albedo.f90:76
ice_albedo::alvdrn
real(kind=dbl_kind), dimension(:,:,:), allocatable, save alvdrn
Definition: ice_albedo.f90:76
ice_therm_vertical::tmlt
real(kind=dbl_kind), dimension(nilyr+1) tmlt
Definition: ice_therm_vertical.f90:78
ice_domain::ilo
integer(kind=int_kind) ilo
Definition: ice_domain.f90:101
ice_constants::p5
real(kind=dbl_kind), parameter p5
Definition: ice_constants.f90:150
ice_constants::puny
real(kind=dbl_kind), parameter puny
Definition: ice_constants.f90:150
ice_constants::tffresh
real(kind=dbl_kind), parameter tffresh
Definition: ice_constants.f90:85
ice_therm_vertical::tridiag_solver
subroutine tridiag_solver(isolve, indxii, indxjj, nmat, sbdiag, diag, spdiag, rhs, xout)
Definition: ice_therm_vertical.f90:2245
ice_calendar::dtice
real(kind=dbl_kind) dtice
Definition: ice_calendar.f90:75
ice_domain::jhi
integer(kind=int_kind) jhi
Definition: ice_domain.f90:101
ice_flux::meltt
real(kind=dbl_kind), dimension(:,:), allocatable, save meltt
Definition: ice_flux.f90:140
ice_state::tsfcn
real(kind=dbl_kind), dimension(:,:,:), allocatable, target, save tsfcn
Definition: ice_state.f90:97
ice_domain
Definition: ice_domain.f90:36
ice_work
Definition: ice_work.f90:38
mod_utils::dbg_sbrio
integer, parameter dbg_sbrio
Definition: mod_utils.f90:70
ice_work::worka
real(kind=dbl_kind), dimension(:,:), allocatable worka
Definition: ice_work.f90:61
ice_constants::p001
real(kind=dbl_kind), parameter p001
Definition: ice_constants.f90:150
ice_atmo
Definition: ice_atmo.f90:36
ice_calendar
Definition: ice_calendar.f90:35
ice_flux::merge_fluxes
subroutine merge_fluxes(ni, strxn, stryn, fsensn, flatn, fswabsn, flwoutn, evapn, Trefn, Qrefn, freshn, fsaltn, fhnetn, fswthrun)
Definition: ice_flux.f90:339
ice_therm_vertical::add_new_snow
subroutine add_new_snow(ni, icells, indxi, indxj, hsn, hsn_new, qsn, Tsf)
Definition: ice_therm_vertical.f90:2952
ice_constants::rhow
real(kind=dbl_kind), parameter rhow
Definition: ice_constants.f90:85
ice_state::vicen
real(kind=dbl_kind), dimension(:,:,:), allocatable, target, save vicen
Definition: ice_state.f90:97
ice_therm_vertical::update_state_vthermo
subroutine update_state_vthermo(ni, icells, indxi, indxj, hin, hsn, qin, qsn, Tsf)
Definition: ice_therm_vertical.f90:3042
ice_therm_vertical::hsnomin
real(kind=dbl_kind), parameter hsnomin
Definition: ice_therm_vertical.f90:73
ice_state::aice
real(kind=dbl_kind), dimension(:,:), allocatable, target, save aice
Definition: ice_state.f90:82
ice_state::eicen
real(kind=dbl_kind), dimension(:,:,:), allocatable, target, save eicen
Definition: ice_state.f90:108
ice_kinds_mod
Definition: ice_kinds_mod.f90:35
ice_atmo::atmo_boundary_layer
subroutine atmo_boundary_layer(ni, sfctype, Tsf, strx, stry, Trf, Qrf, delt, delq)
Definition: ice_atmo.f90:62
ice_model_size::ncat
integer(kind=int_kind), parameter ncat
Definition: ice_model_size.f90:48
ice_constants::c2i
real(kind=dbl_kind), parameter c2i
Definition: ice_constants.f90:150
ice_state::aicen
real(kind=dbl_kind), dimension(:,:,:), allocatable, target, save aicen
Definition: ice_state.f90:97
ice_therm_vertical::thickness_changes
subroutine thickness_changes(ni, icells, indxi, indxj, hin, hlyr, hsn, qin, qsn, mvap, Tbot, fbot, flatn, fhnetn, fsurf, fcondtop, fcondbot, efinal)
Definition: ice_therm_vertical.f90:2349
ice_flux::swvdr
real(kind=dbl_kind), dimension(:,:), allocatable, save swvdr
Definition: ice_flux.f90:91
ice_flux::init_diagnostics
subroutine init_diagnostics
Definition: ice_flux.f90:286
ice_flux::frzmlt
real(kind=dbl_kind), dimension(:,:), allocatable, save frzmlt
Definition: ice_flux.f90:91
ice_calendar::yday
real(kind=dbl_kind) yday
Definition: ice_calendar.f90:75
ice_therm_vertical::ferrmax
real(kind=dbl_kind), parameter ferrmax
Definition: ice_therm_vertical.f90:73
ice_flux::meltb
real(kind=dbl_kind), dimension(:,:), allocatable, save meltb
Definition: ice_flux.f90:140
ice_fileunits
Definition: ice_fileunits.f90:35
ice_constants
Definition: ice_constants.f90:36
ice_flux::shcoef
real(kind=dbl_kind), dimension(:,:), allocatable, save shcoef
Definition: ice_flux.f90:164
ice_flux::strocnyt
real(kind=dbl_kind), dimension(:,:), allocatable, save strocnyt
Definition: ice_flux.f90:55
ice_flux::fsnow
real(kind=dbl_kind), dimension(:,:), allocatable, save fsnow
Definition: ice_flux.f90:91
ice_constants::c1i
real(kind=dbl_kind), parameter c1i
Definition: ice_constants.f90:150
ice_work::workb
real(kind=dbl_kind), dimension(:,:), allocatable workb
Definition: ice_work.f90:61
ice_flux::qa
real(kind=dbl_kind), dimension(:,:), allocatable, save qa
Definition: ice_flux.f90:91
ice_flux::swidf
real(kind=dbl_kind), dimension(:,:), allocatable, save swidf
Definition: ice_flux.f90:91
ice_albedo
Definition: ice_albedo.f90:37
ice_flux::frz_onset
real(kind=dbl_kind), dimension(:,:), allocatable, save frz_onset
Definition: ice_flux.f90:140
ice_constants::ksno
real(kind=dbl_kind), parameter ksno
Definition: ice_constants.f90:118
ice_therm_vertical::ustar_scale
real(kind=dbl_kind) ustar_scale
Definition: ice_therm_vertical.f90:78
ice_therm_vertical::init_thermo_vars
subroutine init_thermo_vars(strxn, stryn, Trefn, Qrefn, fsensn, flatn, fswabsn, flwoutn, evapn, freshn, fsaltn, fhnetn, fswthrun, fsurf, fcondtop, fcondbot, fswint, einit, efinal, mvap)
Definition: ice_therm_vertical.f90:675
ice_itd
Definition: ice_itd.f90:45
ice_therm_vertical::l_brine
logical(kind=log_kind) l_brine
Definition: ice_therm_vertical.f90:91
ice_therm_vertical::rside
real(kind=dbl_kind), dimension(:,:), allocatable, save rside
Definition: ice_therm_vertical.f90:88
ice_therm_vertical::temperature_changes
subroutine temperature_changes(ni, icells, indxi, indxj, hlyr, hsn, qin, Tin, qsn, Tsn, Tsf, Tbot, fbot, fsensn, flatn, fswabsn, flwoutn, fswthrun, fhnetn, fsurf, fcondtop, fcondbot, fswint, einit)
Definition: ice_therm_vertical.f90:1128
ice_therm_vertical::thermo_vertical
subroutine thermo_vertical
Definition: ice_therm_vertical.f90:188
ice_constants::rhoi
real(kind=dbl_kind), parameter rhoi
Definition: ice_constants.f90:85
ice_constants::kice
real(kind=dbl_kind), parameter kice
Definition: ice_constants.f90:118
ice_flux
Definition: ice_flux.f90:36
ice_constants::p333
real(kind=dbl_kind), parameter p333
Definition: ice_constants.f90:150
ice_constants::eps11
real(kind=dbl_kind), parameter eps11
Definition: ice_constants.f90:150
ice_therm_vertical::surface_fluxes
subroutine surface_fluxes(isolve, indxii, indxjj, Tsf, fswsfc, flwoutn, fsensn, flatn, fsurf, dflwout_dT, dfsens_dT, dflat_dT, dfsurf_dT)
Definition: ice_therm_vertical.f90:2146
ice_flux::swidr
real(kind=dbl_kind), dimension(:,:), allocatable, save swidr
Definition: ice_flux.f90:91
ice_flux::mlt_onset
real(kind=dbl_kind), dimension(:,:), allocatable, save mlt_onset
Definition: ice_flux.f90:140
ice_itd::ilyr1
integer(kind=int_kind), dimension(ncat) ilyr1
Definition: ice_itd.f90:62
ice_constants::kappav
real(kind=dbl_kind), parameter kappav
Definition: ice_constants.f90:118
ice_therm_vertical::absorbed_solar
subroutine absorbed_solar(ni, icells, indxi, indxj, hlyr, hsn, fswsfc, fswint, fswthrun, Iabs)
Definition: ice_therm_vertical.f90:1996
ice_therm_vertical::salin
real(kind=dbl_kind), dimension(nilyr+1) salin
Definition: ice_therm_vertical.f90:78
ice_therm_vertical
Definition: ice_therm_vertical.f90:46
ice_albedo::alvdfn
real(kind=dbl_kind), dimension(:,:,:), allocatable, save alvdfn
Definition: ice_albedo.f90:76
ice_therm_vertical::conductivity
subroutine conductivity(icells, indxi, indxj, hlyr, hsn, Tin, Tbot, khi, khs)
Definition: ice_therm_vertical.f90:1887
ice_model_size
Definition: ice_model_size.f90:34
ice_state::vsnon
real(kind=dbl_kind), dimension(:,:,:), allocatable, target, save vsnon
Definition: ice_state.f90:97
ice_constants::cp_ocn
real(kind=dbl_kind), parameter cp_ocn
Definition: ice_constants.f90:85
Generated by   doxygen 1.8.15