forked from NOAA-GFDL/SIS2
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathSIS2_ice_thm.F90
2074 lines (1792 loc) · 95.5 KB
/
SIS2_ice_thm.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
! !
! N-LAYER VERTICAL THERMODYNAMICS !
! !
! References: !
! Hallberg, R., and M. Winton, 2016: The SIS2.0 sea-ice model, in prep. !
! !
! Winton, M., 2011: A conservative non-iterative n-layer sea ice !
! temperature solver, in prep. !
! !
! !
! ->+---------+ <- ts - diagnostic surface temperature ( <= 0C ) !
! / | | !
! hs | snow | <- tsn One snow layer with heat capacity !
! \ | | !
! =>+---------+ !
! / | | !
! / | | <- t1 N salty ice layers with heat capacity !
! / | | !
! / | | <- t2 !
! hi |...ice...| !
! \ | | <- tN-1 !
! \ | | !
! \ | | <- tN !
! \ | | !
! ->+---------+ <- base of ice fixed at seawater freezing temp. !
! !
! Bob Hallberg (Robert.Hallberg@noaa.gov)!
! Mike Winton (Michael.Winton@noaa.gov) !
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
! Enhancement for melt ponds cohabiting ice top with snow layer: !
! !
! ->+--------+ <- ts==0C when pond (net heat here -> pond freeze/melt) !
! / | |--------|<-+ !
! hs | snow | pond | hp - layer has no heat capacity; "surface" energy !
! | | | . balance takes place at ice/snow top !
! \ | | | . !
! =>+--------+--------|<-+ !
! / | | !
! hi | ...ice... | do_pond = true activates scheme !
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
!
! SIS2 "getting-started" melt pond scheme
!
! This melt pond scheme adds a single layer melt pond to each ice thickness
! category. The layer does not have heat capacity. It is assumed to be at
! freshwater freezing temperature and well mixed. All pond surface fluxes
! are communicated directly to its bottom where surface energy balance is
! calculated. The pond layer is advected and redistributed between ice
! thickness categories similarly to the snow and ice layers. As is the case
! with snow, pond cannot exist without an ice layer below. Basic descriptions
! of pond sources and sinks, and pond fraction for radiation treatments follow:
!
! Source water: Surface melting of snow and ice are the source of pond water.
! The runoff scheme keys on the total ice covered area. This is the only
! interaction between category variables in the scheme. The scheme uses r, the
! fraction of melt (r)etained given by the CICE5 scheme (p. 42):
! r = r_min + (r_max-r_min)*ice_area. The controlling namelist
! parameters for r_min and r_max are pond_r_min and pond_r_max.
!
! Freezing sink: The surface energy budget continues to be performed at the
! top of the snow or ice when pond is present but the interface is fixed at
! freezing temperature and the residual between the fluxes on either side of the
! interface is made up with melt or freezing. Freezing in this calculation is a
! sink of pond water.
!
! Freeboard sink: if the pond and snow are sufficiently massive to push the top
! of the sea ice below sea level, pond is dumped into the ocean until the ice
! top is brought back to sea level, or the pond is completely depleted. This
! adjustment follows the CICE5 Hunke "level ice" pond scheme.
!
! Porous-ice sink: No through-ice drainage occurs until the ice average
! temperature exceeds a specified value. The namelist parameter for this is
! pond_porous_temp. Once this limit is exceeded the pond drains to a minimum
! value intended to represent coverage by ponds at sea level. This scheme is
! a placeholder for the mushy-layer thermodynamics to be implemented later.
!
! Pond fraction for radiation: In the snow-free case we assume that pond
! fraction ranges between a specified minimum, where surface cavities below
! sea level are filled with pond water, and a specified maximum value.
! The pond depth and pond fraction are considered to be proportional as was
! found at the SHEBA site and incorporated into the Bailey melt pond
! parameterization. This means that the pond fraction is proportional to
! the square root of the pond volume. The pond volume has a maximum
! determined by the non-negative freeboard requirement. This pond volume
! is associated with the maximum pond fraction, so less pond water is needed
! to cover thinner ice.
!
! New namelist parameters: do_pond, pond_r_max, pond_r_min, pond_porous_temp,
! pond_frac_max, pond_frac_min
!
! New diagnostics: hp, fp, pond_source, pond_sink_freeboard, pond_sink_porous,
! pond_sink_tot <only hp implemented so far; add pond transport diagnostics?>
!
! M. Winton (6/16)
!
module SIS2_ice_thm
use ice_thm_mod, only : get_thermo_coefs
use MOM_EOS, only : EOS_type, EOS_init, EOS_end
use MOM_error_handler, only : SIS_error=>MOM_error, FATAL, WARNING, SIS_mesg=>MOM_mesg
use MOM_file_parser, only : get_param, log_param, read_param, log_version, param_file_type
use MOM_obsolete_params, only : obsolete_logical
implicit none ; private
public :: SIS2_ice_thm_init, SIS2_ice_thm_end, ice_temp_SIS2
public :: ice_resize_SIS2, add_frazil_SIS2, rebalance_ice_layers
public :: get_SIS2_thermo_coefs, ice_thermo_init, ice_thermo_end
public :: Temp_from_Enth_S, Temp_from_En_S, enth_from_TS, enthalpy_from_TS
public :: enthalpy_liquid_freeze, T_Freeze, calculate_T_Freeze, enthalpy_liquid
public :: e_to_melt_TS, energy_melt_enthS
type, public :: ice_thermo_type ; private
real :: Cp_ice ! The heat capacity of ice, in J kg-1 K-1.
real :: Cp_water ! The heat capacity of liquid water in the ice model,
! but not in the brine pockets, in J/(kg K).
real :: Cp_brine ! The heat capacity of liquid water in the brine
! pockets within the ice, in J/(kg K). Cp_brine
! should be set equal to Cp_Water, but for
! algorithmic convenience can be set equal to Cp_ice.
real :: rho_ice, rho_snow, rho_water ! The nominal densities of ice and water in kg m-3.
real :: LI ! The latent heat of fusion, in J kg-1.
real :: Lat_Vapor ! The latent heat of vaporization, in J kg-1.
real :: dTf_dS ! The derivative of the freezing point with salinity,
! in degC per PSU. (dTf_dS is negative.)
real :: enth_liq_0 = 0.0 ! The value of enthalpy for liquid fresh
! water at 0 C, in J kg-1.
real :: enth_unit = 1.0 ! A conversion factor for enthalpy from Joules kg-1.
logical :: slab_ice = .false. ! If true use the very old slab ice thermodynamics,
! with effectively zero heat capacity of ice and snow.
type(EOS_type), pointer :: EOS=>NULL() ! A pointer to the shared MOM6/SIS2
! equation-of-state type. This is here to encourage
! the use of common and consistent thermodynamics
! between the ice and ocean.
end type ice_thermo_type
type, public :: SIS2_ice_thm_CS ; private
! properties of ice, snow, and seawater (NCAR CSM values)
real :: KS ! conductivity of snow, often 0.31 W/(mK)
real :: KI ! conductivity of ice, often 2.03 W/(mK)
real :: temp_ice_freeze ! The freezing temperature of the top ice layer, in C.
real :: temp_range_est ! An estimate of the range of snow and ice temperatures
! that is used to evaluate whether an explicit
! diffusive form of the heat fluxes or an inversion
! based on the layer heat budget is more likely to
! be the most accurate.
real :: h_lo_lim ! hi/hs lower limit for temp. calc.
real :: frazil_temp_offset = 0.5 ! An offset between the temperature with
! which frazil ice forms and the freezing point of
! each sublayer of the ice. This functionality could
! later be accounted for using liq_lim instead.
! In the ice temperature calculation we place a limit to below (salinity
! dependent) freezing point on the prognosed temperatures. For ice_resize
! it is better to make a slightly more restrictive limit that requires the
! temperature to be such that the brine content is less than "liq_lim" of
! the total mass. That is T_f/T < liq_lim implying T<T_f/liq_lim
real :: liq_lim = .99
logical :: do_pond = .false. ! activate melt pond scheme - mw/new
! mw/new - these melt pond control data are temporarily placed here
real :: tdrain = -0.8 ! if average ice temp. > tdrain, drain pond
real :: r_min_pond = 0.15 ! pond retention of meltwater
real :: r_max_pond = 0.9 ! see CICE5 doc
real :: max_pond_frac = 0.5 ! pond water beyond this is dumped
real :: min_pond_frac = 0.2 ! ponds below sea level don't drain
! mw/new - end of melt pond control data
end type SIS2_ice_thm_CS
contains
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
!> SIS2_ice_thm_init initializes the control structure for the ice thermodynamic
!! update code.
subroutine SIS2_ice_thm_init(param_file, CS)
type(param_file_type), intent(in) :: param_file
type(SIS2_ice_thm_CS), pointer :: CS
! This include declares and sets the variable "version".
#include "version_variable.h"
character(len=40) :: mod = "SIS2_ice_thm (updates)" ! This module's name.
if (.not.associated(CS)) allocate(CS)
call log_version(param_file, mod, version, &
"This sub-module does updates of the sea-ice due to thermodynamic changes.")
call get_param(param_file, mod, "SNOW_CONDUCTIVITY", CS%Ks, &
"The conductivity of heat in snow.", units="W m-1 K-1", &
default=0.31)
call get_param(param_file, mod, "ICE_CONDUCTIVITY", CS%Ki, &
"The conductivity of heat in ice.", units="W m-1 K-1", &
default=2.03)
call get_param(param_file, mod, "MIN_H_FOR_TEMP_CALC", CS%h_lo_lim, &
"The minimum ice thickness at which to do temperature \n"//&
"calculations.", units="m", default=0.0)
call get_param(param_file, mod, "DO_POND", CS%do_pond, &
"If true, calculate melt ponds and use them for\n"//&
"shortwave radiation calculation.", default=.false.)
call get_param(param_file, mod, "TDRAIN", CS%tdrain, &
"Melt ponds drain to sea level when ice average temp.\n"//&
"exceeds TDRAIN (stand-in for mushy layer thermo)", default=-0.8)
call get_param(param_file, mod, "R_MIN_POND", CS%r_min_pond, &
"Minimum retention rate of surface water sources in melt pond\n"//&
"(retention scales linearly with ice cover)", default=0.15)
call get_param(param_file, mod, "R_MAX_POND", CS%r_max_pond, &
"Maximum retention rate of surface water sources in melt pond\n"//&
"(retention scales linearly with ice cover)", default=0.9)
call get_param(param_file, mod, "MIN_POND_FRAC", CS%min_pond_frac, &
"Minimum melt pond cover (by ponds at sea level)\n"//&
"pond drains to this when ice is porous.", default=0.2)
call get_param(param_file, mod, "MAX_POND_FRAC", CS%max_pond_frac, &
"Maximum melt pond cover - associated with pond volume\n"//&
"that suppresses ice top to waterline", default=0.5)
call get_param(param_file, mod, "ICE_TEMP_RANGE_ESTIMATE", CS%temp_range_est,&
"An estimate of the range of snow and ice temperatures \n"//&
"that is used to evaluate whether an explicit diffusive \n"//&
"form of the heat fluxes or an inversion based on the \n"//&
"layer heat budget is more likely to be more accurate. \n"//&
"Setting this to 0 causes the explicit diffusive form. \n"//&
"to always be used.", units="degC", default=40.0)
CS%temp_range_est = abs(CS%temp_range_est)
call obsolete_logical(param_file, "OLD_ICE_HEAT_CAPACITY", warning_val=.false.)
end subroutine SIS2_ice_thm_init
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
! ice_temp_SIS2 - A subroutine that calculates the snow and ice enthalpy !
! changes due to surface forcing. !
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
! m_pond - mw/new
subroutine ice_temp_SIS2(m_pond, m_snow, m_ice, enthalpy, sice, sh_T0, B, sol, tfw, fb, &
tsurf, dtt, NkIce, tmelt, bmelt, CS, ITV, check_conserve)
real, intent(in ) :: m_pond ! pond mass per unit area (kg m-2)
real, intent(in ) :: m_snow ! snow mass per unit area (H, usually kg m-2)
real, intent(in ) :: m_ice ! ice mass per unit area (H, usually kg m-2)
real, dimension(0:NkIce) , &
intent(inout) :: enthalpy ! The enthalpy of each layer in a column of
! snow and ice, in enth_unit (J kg-1).
real, dimension(NkIce), &
intent(in) :: Sice ! ice salinity by layer (g/kg)
real, intent(in ) :: sh_T0 ! net surface heat flux (+ up) at ts=0 (W/m^2)
real, intent(in ) :: B ! d(sfc heat flux)/d(ts) [W/(m^2 deg-C)]
real, dimension(0:NkIce), &
intent(in) :: sol ! Solar heating of the snow and ice layers (W m-2)
real, intent(in ) :: tfw ! seawater freezing temperature (deg-C)
real, intent(in ) :: fb ! heat flux upward from ocean to ice bottom (W/m^2)
real, intent( out) :: tsurf ! surface temperature (deg-C)
real, intent(in ) :: dtt ! timestep (sec)
integer, intent(in ) :: NkIce ! The number of ice layers.
real, intent(inout) :: tmelt ! accumulated top melting energy (J/m^2)
real, intent(inout) :: bmelt ! accumulated bottom melting energy (J/m^2)
type(SIS2_ice_thm_CS), intent(in) :: CS
type(ice_thermo_type), intent(in) :: ITV ! The ice thermodynamic parameter structure.
logical, optional, intent(in) :: check_conserve ! If true, check for local heat conservation.
!
! variables for temperature calculation [see Winton (1999) section II.A.]
! note: here equations are multiplied by hi to improve thin ice accuracy
!
real :: A ! Net downward surface heat flux from the atmosphere at 0C (W/m^2)
! real, dimension(0:NkIce) :: &
! temp_est, & ! An estimated snow and ice temperature, in degC.
! temp_IC, & ! The temperatures of the snow and ice based on the initial
! ! enthalpy, in degC.
! temp_new ! The updated temperatures, in degC.
real, dimension(0:NkIce) :: temp_est ! An estimated snow and ice temperature, in degC.
real, dimension(0:NkIce) :: temp_IC ! The temperatures of the snow and ice based on the initial
! enthalpy, in degC.
real, dimension(0:NkIce) :: temp_new ! The updated temperatures, in degC.
real, dimension(NkIce) :: tfi ! The ice freezing temperatures, in degC.
real :: mL_ice ! The mass-per-unit-area of each ice layer in kg m-2 (not H).
real :: mL_snow ! The mass-per-unit-area of each snow layer in kg m-2 (not H).
real :: e_extra
real, dimension(0:NkIce) :: m_lay ! Masses of all layers in kg m-2.
real :: enth_fp ! The enthalpy at the freezing point (solid for fresh ice).
real :: kk, k10, k0a, k0skin
real :: I_bb, b_denom_1
real :: comp_rat ! The complement of rat, going from 0 to 1.
real :: tsf ! The surface freezing temperature in degC.
real :: k0a_x_ta, tsno_est, salt_part, rat
real :: tsurf_est ! An estimate of the surface temperature in degC.
real, dimension(0:NkIce+1) :: cc ! Interfacial coupling coefficients.
real, dimension(0:NkIce) :: bb ! Effective layer heat capacities.
real, dimension(0:NkIce) :: cc_bb ! Remaining coupling ratios.
real, dimension(-1:NkIce) :: heat_flux_int ! The downward heat fluxes at the
! interfaces between layers, in W m-2.
! heat_flux_int uses the index convention from
! MOM6 that interface K is below layer k.
real :: I_liq_lim ! The inverse of CS%liq_lim.
real :: heat_flux_err_rat
real :: col_enth1, col_enth2, col_enth2b, col_enth3
real :: d_e_extra, e_extra_sum
real :: tflux_bot, tflux_bot_diff, tflux_bot_resid
real :: tfb_diff_err, tfb_resid_err
real :: tflux_sfc, sum_sol, d_tflux_bot
real :: hsnow_eff
real :: snow_temp_new, snow_temp_max
real :: hL_ice_eff
real :: enth_liq_lim
real :: enth_prev
real :: rho_ice ! The nominal density of sea ice in kg m-3.
real :: rho_snow ! The nominal density of snow in kg m-3.
real :: Cp_ice ! The heat capacity of ice, in J kg-1 K-1.
real :: Cp_brine ! The heat capacity of liquid water in the brine pockets,
! in J kg-1 K-1.
real :: Lat_fus ! The latent heat of fusion, in J kg-1.
real :: enth_unit ! A conversion factor for enthalpy from Joules kg-1.
real :: I_enth_unit ! The inverse of enth_unit.
logical :: col_check
integer :: k
col_check = .false. ; if (present(check_conserve)) col_check = check_conserve
temp_IC(0) = temp_from_En_S(enthalpy(0), 0.0, ITV)
call temp_from_Enth_S(enthalpy(1:), sice(1:), temp_IC(1:), ITV)
call get_SIS2_thermo_coefs(ITV, enthalpy_units=enth_unit, rho_ice=rho_ice, rho_snow=rho_snow, &
Cp_Ice=Cp_ice, Latent_Fusion=Lat_fus, Cp_Brine=Cp_Brine)
A = -sh_T0
I_enth_unit = 1.0 / enth_unit
mL_ice = m_ice / NkIce ! ice mass per unit area of each layer
mL_snow = m_snow ! snow mass per unit area (in kg m-2).
call calculate_T_Freeze(sice, tfi, ITV) ! freezing temperature of ice layers
! Set the effective thickness of each ice and snow layer, limited to avoid
! instabilities for thin layers.
hL_ice_eff = max(mL_ice / rho_ice, CS%H_LO_LIM)
hsnow_eff = mL_snow / rho_snow + max(1e-35, 1e-20*CS%H_LO_LIM)
kk = CS%KI/hL_ice_eff ! full ice layer conductivity
tsf = tfi(1) ! surface freezing temperature
if (mL_snow>0.0) tsf = 0.0
k10 = 2.0*(CS%KS*CS%KI) / (hL_ice_eff*CS%KS + hsnow_eff*CS%KI) ! coupling ice layer 1 to snow
k0a = (CS%KS*B) / (0.5*B*hsnow_eff + CS%KS) ! coupling snow to "air"
k0skin = 2.0*CS%KS / hsnow_eff
k0a_x_ta = (CS%KS*A) / (0.5*B*hsnow_eff + CS%KS) ! coupling times "air" temperture
enth_liq_lim = Enth_from_TS(0.0, 0.0, ITV)
! Determine the enthalpy for conservation checks.
m_lay(0) = mL_snow ; do k=1,NkIce ; m_lay(k) = mL_ice ; enddo
if (col_check) then
col_enth1 = 0.0
do k=0,NkIce ; col_enth1 = col_enth1 + m_lay(k)*enthalpy(k) ; enddo
endif
e_extra_sum = 0.0
! First get a non-conservative estimate with a linearized fully implicit
! treatment of layer coupling.
! Determine the effective layer heat capacities.
! bb = dheat/dTemp, as derived from a linearization of the enthalpy equation
! but using the value for ice near the melt point if temp_IC is above T_fr.
bb(0) = mL_snow* Cp_ice
do k=1,NkIce ! Store the heat capacity term in bb.
if (tfi(k) >= 0.0) then ! This is pure ice.
bb(k) = mL_ice * Cp_ice
elseif (temp_IC(k) < tfi(k)) then
bb(k) = mL_ice * (Cp_ice - (tfi(k) / temp_IC(k)**2) * &
(Lat_Fus - (Cp_brine-Cp_ice) * temp_IC(k)) )
! Or more generally:
! S_Sf = sice(k) / calculate_S_freeze(temp_IC(k))
! bb(k) = mL_ice * (Cp_ice + dSf_dT*Lat_fus + S_Sf * &
! ((Cp_brine-Cp_ice))
else ! Use the value when temp_IC = tfi.
bb(k) = mL_ice * (Cp_brine - Lat_fus / tfi(k))
endif
enddo
! The following expressions could be modified to permit there to be
! variable diffusivities in the ice/brine mixtures.
cc(0) = k0a*dtt ! Atmosphere-snow coupling
cc(1) = k10*dtt ! Snow-ice coupling
do k=2,NkIce ; cc(k) = kk*dtt ; enddo
cc(NkIce+1) = 2*kk*dtt ! bottom of ice coupling with ocean at freezing point.
! This is a version of a tridiagonal solver where the relationship between
! the diagonal elements and the coupling between layers has been used to
! ensure that there are no differences ever taken, so there will be minimal
! truncation errors. This closely follows schemes that have previously been
! found to work very well in GOLD and MOM6.
! Go UP the ice column.
b_denom_1 = (bb(NkIce) + cc(NkIce+1))
I_bb = 1.0 / (b_denom_1 + cc(NkIce))
temp_est(NkIce) = ( (sol(NkIce)*dtt + bb(NkIce)*temp_IC(NkIce)) + &
cc(NkIce+1)*tfw ) * I_bb
comp_rat = b_denom_1 * I_bb
cc_bb(NkIce) = cc(NkIce) * I_bb
do k=NkIce-1,1,-1
b_denom_1 = bb(k) + comp_rat*cc(k+1)
I_bb = 1.0 / (b_denom_1 + cc(k))
temp_est(k) = ((sol(k)*dtt + bb(k)*temp_IC(k)) + cc(k+1)*temp_est(k+1)) * I_bb
comp_rat = b_denom_1 * I_bb ! 1.0 >= comp_rat >= 0.0
cc_bb(k) = cc(k) * I_bb
end do
b_denom_1 = bb(0) + comp_rat*cc(1)
I_bb = 1.0 / (b_denom_1 + cc(0))
! This is a complete calculation of temp_est(0), assuming that the surface
! flux is given by A - B*tsurf, with tsurf estimated as part of this calculation.
temp_est(0) = (((sol(0)*dtt + bb(0)*temp_IC(0)) + k0a_x_ta*dtt) + cc(1)*temp_est(1)) * I_bb
! Diagnose the surface skin temperature by matching the diffusive fluxes in
! the snow with the atmospheric fluxes. I.e. solve the following for tsurf_est:
! (A - B*tsurf_est) = k0skin * (tsurf_est - temp_est(0))
tsurf_est = (A + k0skin*temp_est(0)) / (B + k0skin)
if (tsurf_est > tsf .or. m_pond > 0.0 ) then ! mw/new - liq. h2o @ sfc
! also pins temp at freezing
! The surface is melting, set tsurf to melt temp. and recalculate I_bb.
tsurf_est = tsf
! cc(0) = k0skin*dtt
I_bb = 1.0 / (b_denom_1 + k0skin*dtt)
temp_est(0) = min(tsf, &
(((sol(0)*dtt + bb(0)*temp_IC(0)) + k0skin*dtt*tsf) + cc(1)*temp_est(1)) * I_bb)
endif
! Go back DOWN the ice column to get the estimated temperatures, subject to the
! limitation that all temperatures are assumed to be at or below freezing.
do k=1,NkIce
temp_est(k) = min(temp_est(k) + cc_bb(k) * temp_est(k-1), tfi(k))
end do
!
! Quasi-conservative iterative pass going UP the ice column
!
temp_est(NkIce) = laytemp_SIS2(mL_ice, tfi(NkIce), sol(NkIce) + kk*(2*tfw+temp_est(NkIce-1)), &
3*kk, temp_IC(NkIce), enthalpy(NkIce), sice(NkIce), dtt, ITV)
do k=NkIce-1,2,-1
temp_est(k) = laytemp_SIS2(mL_ice, tfi(k), sol(k) + kk*(temp_est(k-1)+temp_est(k+1)), &
2*kk, temp_IC(k), enthalpy(k), sice(k), dtt, ITV)
enddo
temp_est(1) = laytemp_SIS2(mL_ice, tfi(1), sol(1) + (kk*temp_est(2) + k10*temp_est(0)), &
kk + k10, temp_IC(1), enthalpy(1), sice(1), dtt, ITV)
! Calculate the bulk snow temperature and surface skin temperature together.
temp_est(0) = laytemp_SIS2(mL_snow, 0.0, sol(0) + (k10*temp_est(1)+k0a_x_ta), &
k10+k0a, temp_IC(0), enthalpy(0), 0.0, dtt, ITV)
tsurf = (A + k0skin*temp_est(0)) / (B + k0skin) ! diagnose surface skin temp.
!
! The following conservative update pass going DOWN the ice column is where
! the layer enthalpies are actually updated and extra heat that should drive
! melting is accumulated and stored.
!
! Pre-calculate a conversion factor that can be used to figure out whether it
! is more accurate to use the explicit form for the fluxes or to invert enthalpy
! conservation for the fluxes, depending on the relative layer masses and the
! strength of the coupling between layers. The expression below is a simplified
! form that assume that temperatures are measured in Celsius and are less than
! about CS%temp_err_est, and that the heat budget has 4 terms of comparable
! magnitude to the enthalpy content of the layer. If there were a larger
! tolerance for the iterative estimate of a new temperature, that would go into
! the numerator but not the denominator. The present form is based on the
! assumption that floating-point roundoff dominates the errors.
heat_flux_err_rat = 0.7071 * dtt * (CS%temp_range_est) / &
(CS%temp_range_est * Cp_ice + Lat_fus)
e_extra = 0.0
if (tsurf > tsf .or. m_pond > 0.0) then
! The surface is melting: update enthalpy with the surface at melt temp.
! mw/new - liq h2o at surface also pins temp at freezing
tsurf = tsf
! Accumulate surface melt energy.
if (mL_snow>0.0) then
heat_flux_int(-1) = k0skin * tsf
heat_flux_int(0) = -k10*temp_est(1)
call update_lay_enth(mL_snow, 0.0, enthalpy(0), heat_flux_int(-1), &
sol(0), heat_flux_int(0), -k0skin, k10, dtt, &
heat_flux_err_rat, ITV, e_extra)
tmelt = tmelt + e_extra + dtt*((A-B*tsf) - heat_flux_int(-1))
tflux_sfc = dtt*heat_flux_int(-1)
e_extra_sum = e_extra_sum + e_extra
else
! There is no snow mass, so just convert tsnow = tsf to enthalpy.
enthalpy(0) = enth_from_TS(tsf, 0.0, ITV)
heat_flux_int(0) = k10*(tsf - temp_est(1))
heat_flux_int(-1) = heat_flux_int(0)
! Replace use tsurf in the calculation of tmelt
tmelt = tmelt + dtt*((sol(0)+(A-B*tsf)) - heat_flux_int(0))
tflux_sfc = dtt*heat_flux_int(0)
endif
else
heat_flux_int(-1) = k0a_x_ta
heat_flux_int(0) = -k10*temp_est(1)
snow_temp_max = (tsf*(B + k0skin) - A) / k0skin
call update_lay_enth(mL_snow, 0.0, enthalpy(0), heat_flux_int(-1), &
sol(0), heat_flux_int(0), -k0a, k10, dtt, &
heat_flux_err_rat, ITV, e_extra, &
temp_new=snow_temp_new, temp_max=snow_temp_max)
tsurf = (A + k0skin*snow_temp_new) / (B + k0skin) ! diagnose surface skin temp.
! This is equivalent to, but safer than, tsurf = (A - heat_flux_int(-1)) / B
tflux_sfc = dtt*heat_flux_int(-1)
e_extra_sum = e_extra_sum + e_extra
tmelt = tmelt + e_extra
endif
do k=1,NkIce-1 ! The flux from above is fixed, only have downward feedback
heat_flux_int(K) = -kk*temp_est(k+1)
call update_lay_enth(mL_ice, Sice(k), enthalpy(k), heat_flux_int(K-1), &
sol(k), heat_flux_int(K), 0.0, kk, dtt, &
heat_flux_err_rat, ITV, e_extra)
e_extra_sum = e_extra_sum + e_extra
if (k <= NkIce/2) then ; tmelt = tmelt + e_extra
else ; bmelt = bmelt + e_extra ; endif
enddo
heat_flux_int(NkIce) = -2.0*kk*tfw
call update_lay_enth(mL_ice, Sice(NkIce), enthalpy(NkIce), heat_flux_int(NkIce-1), &
sol(NkIce), heat_flux_int(NkIce), 0.0, 2.0*kk, dtt, &
heat_flux_err_rat, ITV, e_extra)
e_extra_sum = e_extra_sum + e_extra
bmelt = bmelt + e_extra
!
! END of the conservative update of enthalpy.
!
if (col_check) then
col_enth2 = e_extra_sum*enth_unit ; sum_sol = 0.0 ; col_enth2b = 0.0
do k=0,NkIce
col_enth2 = col_enth2 + m_lay(k)*enthalpy(k)
col_enth2b = col_enth2b + m_lay(k)*enthalpy(k)
sum_sol = sum_sol + dtt*sol(k)
enddo
! tflux_bot_resid and tflux_bot_diff are two mathematically equivalent
! estimates of the heat flux at the base of the ice.
tflux_bot_resid = (col_enth2 - col_enth1) - (sum_sol + tflux_sfc)
tflux_bot_diff = -heat_flux_int(NkIce)*dtt
! Estimate the errors with these two expressions from 64-bit roundoff.
tfb_diff_err = 1e-15*2.0*kk*dtt * sqrt(tfw**2 + 10.0**2) ! The -10 deg is arbitrary but good enough?
tfb_resid_err = 1e-15*sqrt(col_enth2**2 + col_enth1**2 + sum_sol**2 + tflux_sfc**2)
d_tflux_bot = tflux_bot_diff - tflux_bot_resid
if (abs(d_tflux_bot) > 1.0e-9*(abs(tflux_bot_resid) + abs(tflux_bot_diff) + &
abs(col_enth2 - col_enth1))) then
d_tflux_bot = tflux_bot_diff - tflux_bot_resid
endif
endif
tflux_bot = -heat_flux_int(NkIce)*dtt
! Accumulate the difference between the ocean's heat flux to the ice-ocean
! interface and the sea-ice heat flux to evaluate bottom melting/freezing.
bmelt = bmelt + (dtt*fb - tflux_bot)
e_extra_sum = 0.0
if (enthalpy(0) > enth_liq_lim) then ! put excess snow energy into top melt.
e_extra = (enthalpy(0) - enth_liq_lim) * mL_snow * I_enth_unit
tmelt = tmelt + e_extra
e_extra_sum = e_extra_sum + e_extra
enthalpy(0) = enth_liq_lim
endif
I_liq_lim = 1.0 / CS%liq_lim
do k=1,NkIce
enth_liq_lim = enth_from_TS(tfi(k)*I_liq_lim, sice(k), ITV)
if (enthalpy(k) > enth_liq_lim) then
! Put excess energy into the closer of top or bottom melt.
e_extra = (enthalpy(k) - enth_liq_lim) * mL_ice * I_enth_unit
e_extra_sum = e_extra_sum + e_extra
enthalpy(k) = enth_liq_lim
if (k<=NkIce/2) then
tmelt = tmelt+e_extra
else
bmelt = bmelt+e_extra
endif
endif
enddo
if (col_check) then
col_enth3 = 0.0
do k=0,NkIce ; col_enth3 = col_enth3 + m_lay(k)*enthalpy(k) ; enddo
d_e_extra = col_enth3 - (col_enth2b - e_extra_sum*enth_unit)
if (abs(d_e_extra) > 1.0e-12*(abs(col_enth3) + abs(col_enth2b) + &
abs(e_extra_sum*enth_unit))) then
d_e_extra = (col_enth3 - col_enth2b) - e_extra_sum*enth_unit
endif
endif
call ice_check(mL_snow, NkIce*mL_ice, enthalpy, sice, NkIce, &
"at end of ice_temp_SIS2", ITV, bmelt=bmelt, tmelt=tmelt, t_sfc=tsurf)
end subroutine ice_temp_SIS2
!
! laytemp_SIS2 - implicit calculation of new layer temperature
!
function laytemp_SIS2(m, T_fr, f, b, tp, enth, salin, dtt, ITV) result (new_temp)
real :: new_temp
real, intent(in) :: m ! mass of ice - kg/m2
real, intent(in) :: T_fr ! ice freezing temp. (determined by salinity)
real, intent(in) :: f ! Inward forcing - W/m2
real, intent(in) :: b ! response of outward heat flux to local temperature - W/m2/K
real, intent(in) :: tp ! prior step temperature
real, intent(in) :: enth ! prior step enthalpy
real, intent(in) :: salin ! ice salinity
real, intent(in) :: dtt ! timestep in s.
type(ice_thermo_type), intent(in) :: ITV ! The ice thermodynamic parameter structure.
real :: T_g ! The latest best guess at Temp, in deg C.
real :: T_deriv ! The value of Temp at which to evaluate dErr_dT, in deg C.
real :: T_max, T_min ! Bracketing temperatures, in deg C.
real :: Err ! The enthalpy at T_guess, in J kg-1.
real :: Err_Tmin, Err_Tmax ! The errors at T_max and T_min, in J m-2.
real :: T_prev ! The previous value of T_g, in deg C.
real :: dErr_dT ! The partial derivative of Err with T_g, in J m-2 C-1.
real :: Enth_tol = 1.0e-15 ! The fractional Enthalpy difference tolerance for convergence.
real :: TfmxdCp_BI
real :: E0 ! Starting heat relative to salinity dependent freezing.
real :: AA, BB, CC
real :: Cp_ice ! The heat capacity of ice, in J kg-1 K-1.
real :: Cp_brine ! The heat capacity of liquid water in the brine pockets,
! in J kg-1 K-1.
real :: Cp_water ! The heat capacity of liquid water in the ice model,
! but not in the brine pockets, in J kg-1 K-1.
real :: LI ! The latent heat of fusion, in J kg-1.
integer :: itt
! real :: T_itt(20), dTemp(20), Err_itt(20)
call get_SIS2_thermo_coefs(ITV, Cp_Ice=Cp_ice, Latent_Fusion=LI, &
Cp_Brine=Cp_Brine, Cp_Water=Cp_Water)
if ( T_fr == 0.0 ) then
! For fresh water, avoid the degeneracy of the enthalpy-temperature
! relationship by extending the linear expression for frozen water, then
! limit it later to be at or below freezing.
!
! m * {Cp_Ice} * (tn-tp) = dtt * (f - b*tn)
!
new_temp = (m*Cp_Ice*tp + f*dtt) / (m*Cp_Ice + b*dtt) ! = -BB/AA
else
if (tp >= T_fr) then
E0 = Cp_Water*(tp - T_fr) ! >= 0
else
E0 = Cp_Ice*(tp - T_fr) - LI*(1 - T_fr/tp) ! < 0
if (Cp_ice /= Cp_brine) E0 = E0 + (Cp_brine - Cp_ice) * T_fr*log(tp/T_fr)
endif
! Determine whether the new solution will be above or below freezing.
if (m*E0 + dtt * (f - b*T_fr) >= 0) then
! This layer will be completely melted, so return the freezing value.
! new_temp = T_fr + (m*E0 + dtt* (f - b*T_fr)) / (Cp_water*m + dtt*b)
new_temp = T_fr
else
! This layer will be partly melted.
! Solve a quadratic equation for the new layer temperature, tn:
!
! m * {Cp_Ice-LI*T_fr/(tn*tp)} * (tn-tp) = dtt * (f - b*tn)
! m * {{Cp_Ice*(tn - T_fr) - LI*(1 - T_fr/tn)} - E0} = dtt * (f - b*tn)
! En0(tn) = Cp_Ice*(tn - T_fr) - LI*(1 - T_fr/tn)
!
AA = m*Cp_Ice + b*dtt
BB = -(m*((E0 + LI) + Cp_Ice*T_fr) + f*dtt)
CC = m*LI*T_fr
! This form avoids round-off errors.
if (BB >= 0) then
new_temp = -(BB + sqrt(BB*BB - 4*AA*CC)) / (2*AA)
else
new_temp = (2*CC) / (-BB + sqrt(BB*BB - 4*AA*CC))
endif
if (Cp_ice /= Cp_brine) then
! At this point, new_temp is just a good starting guess that needs to be iterated to convergence.
! Solve the following expression for the new layer temperature, tn:
!
TfmxdCp_BI = T_fr*m*(Cp_Brine-Cp_Ice)
! Err = m * ((Cp_Ice*(tn - T_fr) - LI*(1 - T_fr/tn)) - &
! TfmxdCp_BI*(log(T_fr/tn)) - E0) - dtt * (f - b*tn)
! Err = -(m*((E0 + LI) + Cp_Ice*T_fr) + f*dtt) + &
! m * (Cp_Ice*tn + LI*(T_fr/tn)) - TfmxdCp_BI*(log(T_fr/tn))) + dtt*b*tn
! This might be a good enough first guess that bracketing is unnecessary
! but it is better to play it safe...
T_g = new_temp
Err = BB + ((m * (Cp_Ice*t_g + LI*(T_fr/t_g)) - &
TfmxdCp_BI*(log(T_fr/t_g))) + dtt*b*t_g)
if (Err <= 0.0) then
T_min = T_g ; Err_Tmin = Err
T_max = T_fr ; Err_Tmax = BB + ((m * (Cp_Ice*T_fr + LI)) + dtt*b*T_fr)
else
T_max = T_g ; Err_Tmax = Err
T_min = -273.15
Err_Tmin = BB + ((m * (Cp_Ice*T_min + LI*(T_fr/T_min)) - &
TfmxdCp_BI*(log(T_fr/T_min))) + dtt*b*T_min)
endif
! T_itt(:) = 0.0 ; dTemp(:) = 0.0 ; Err_itt(:) = 0.0
do itt=1,20 ! Note that 3 or 4 iterations usually are enough.
Err = BB + ((m * (Cp_Ice*t_g + LI*(T_fr/t_g)) - &
TfmxdCp_BI*(log(T_fr/t_g))) + dtt*b*t_g)
! T_itt(itt) = T_g ; Err_itt(itt) = Err
if (abs(Err) <= Enth_tol*(abs(BB) + m*LI + abs(dtt*b*T_fr))) then
new_temp = T_g ; exit
elseif (Err < 0.0) then
T_min = T_g ; Err_Tmin = Err
else
T_max = T_g ; Err_Tmax = Err
endif
! Use the more efficient Newton's method of McDougall & Witherspoon (2014),
! Appl. Math. Lett., 29, 20-25.
T_deriv = T_g
if (itt > 1) then ! Reuse the estimate of dT_dEn from the last iteration.
if ((dErr_dT*T_g - 0.5*Err > dErr_dT*T_min) .and. &
(dErr_dT*T_g - 0.5*Err < dErr_dT*T_max)) &
T_deriv = T_g - 0.5*Err / dErr_dT
endif
dErr_dt = m * (Cp_Ice - LI*T_fr/(t_deriv**2)) + (TfmxdCp_BI / t_deriv + b*dtt) ! >= 0.0
T_prev = T_g
if ((dErr_dT*T_g - Err > dErr_dT*T_min) .and. &
(dErr_dT*T_g - Err < dErr_dT*T_max)) then
T_g = T_g - Err / dErr_dT
else
T_g = (Err_Tmax * T_min - Err_Tmin * T_max) / &
(Err_Tmax - Err_Tmin)
endif
! dTemp(itt) = T_g - T_prev
enddo
new_temp = T_g ! Use the best guess.
! write (*,'("T_itt = ",8F14.8)') T_itt(1:8)
! write (*,'("Err_itt = ",8(1PE14.6))') Err_itt(1:8)
! write (*,'("dTemp = ",8(1Pe12.4))') dTemp(1:8)
endif
endif
endif
! Only return temperatures that are at or below the freezing point.
new_temp = min(new_temp, T_fr)
end function laytemp_SIS2
!
! update_lay_enth - implicit calculation of new layer enthalpy
!
subroutine update_lay_enth(m_lay, sice, enth, ftop, ht_body, fbot, dftop_dT, &
dfbot_dT, dtt, hf_err_rat, ITV, extra_heat, temp_new, temp_max)
real, intent(in) :: m_lay ! This layers mass of ice in kg/m2
real, intent(in) :: sice ! ice salinity in g/kg
real, intent(inout) :: enth ! ice enthalpy in enth_units (proportional to J kg-1).
real, intent(inout) :: ftop ! Downward heat flux atop the layer in W/m2 at T = 0 C, or
! the prescribed heat flux if dftop_dT = 0.
real, intent(in) :: ht_body ! Body forcing to layer in W/m2
real, intent(inout) :: fbot ! Downward heat below the layer in W/m2 at T = 0 C.
real, intent(in) :: dftop_dT ! The linearization of ftop with layer temperature in W m-2 K-1.
real, intent(in) :: dfbot_dT ! The linearization of fbot with layer temperature in W m-2 K-1.
real, intent(in) :: dtt ! The timestep in s.
real, intent(in) :: hf_err_rat ! A conversion factor for comparing the errors
! in explicit and implicit estimates of the updated
! heat fluxes, in (kg m-2) / (W m-2 K-1).
type(ice_thermo_type), intent(in) :: ITV ! The ice thermodynamic parameter structure.
real, intent(out) :: extra_heat ! The heat above the melt point, in J.
real, optional, intent(out) :: temp_new ! The new temperature, in degC.
real, optional, intent(in) :: temp_max ! The maximum new temperature, in degC.
real :: htg ! The rate of heating of the layer in W m-2.
real :: new_temp ! The new layer temperature, in degC.
real :: max_temp ! The maximum new layer temperature, in degC.
real :: max_enth ! The maximum new layer enthalpy, in degC.
real :: fb ! The negative of the dependence of layer heating on
! temperature, in W m-2 K-1. fb > 0.
real :: extra_enth ! Excess enthalpy above the melt point, in kg enth_units.
real :: enth_in ! The initial enthalpy, in enth_units.
real :: enth_fp ! The enthalpy at the freezing point, in enth_units.
real :: AA, BB, CC ! Temporary variables used to solve a quadratic equation.
real :: dtEU ! The timestep times the unit conversion from J to Enth_units, in s?
real :: dT_dEnth ! The partial derivative of temperature with enthalpy,
! in units of K / Enth_unit.
real :: En_J ! The enthalpy in Joules with 0 offset for liquid at 0 C.
real :: T_fr ! Ice freezing temperature (determined by bulk salinity) in deg C.
real :: fbot_in, ftop_in ! Input values of fbot and ftop in W m-2.
real :: dflux_dtot_dT ! A temporary work array in units of degC.
real :: T_g ! The latest best guess at Temp, in deg C.
real :: T_deriv ! The value of Temp at which to evaluate dErr_dT, in deg C.
real :: T_max, T_min ! Bracketing temperatures, in deg C.
real :: Err ! The enthalpy at T_guess, in J kg-1.
real :: Err_Tmin, Err_Tmax ! The errors at T_max and T_min, in J m-2.
real :: T_prev ! The previous value of T_g, in deg C.
real :: dErr_dT ! The partial derivative of Err with T_g, in J m-2 C-1.
real :: Enth_tol = 1.0e-15 ! The fractional Enthalpy difference tolerance for convergence.
real :: TfxdCp_WI, TfxdCp_BI, Err_Tind
real :: Cp_ice ! The heat capacity of ice, in J kg-1 K-1.
real :: Cp_brine ! The heat capacity of liquid water in the brine pockets,
! in J kg-1 K-1.
real :: Cp_water ! The heat capacity of liquid water in the ice model,
! but not in the brine pockets, in J kg-1 K-1.
real :: LI ! The latent heat of fusion, in J kg-1.
real :: enth_unit ! A conversion factor for enthalpy from Joules kg-1.
! real :: Enth_liq_0 ! The enthalpy of liquid water at 0C.
integer :: itt
! real :: T_itt(20), dTemp(20), Err_itt(20)
call get_SIS2_thermo_coefs(ITV, enthalpy_units=enth_unit, Latent_Fusion=LI, &
Cp_Ice=Cp_ice, Cp_Brine=Cp_Brine, Cp_Water=Cp_Water)
! Solve m_lay*(enth - enth_in) + extra_heat = dt * (ht_body + ftop - fbot)
! ftop = ftop_in + temp*dftop_dT
! fbot = fbot_in + temp*dfbot_dT
! subject to enth <= enth_fp and extra_heat >= 0
ftop_in = ftop ; fbot_in = fbot
htg = (ht_body + ftop_in) - fbot_in
fb = -(dftop_dT - dfbot_dT) ! = -dhgt_dt > 0
extra_heat = 0.0 ; extra_enth = 0.0
if (sice > 0.0) then
T_fr = T_freeze(sice, ITV)
enth_fp = enthalpy_liquid_freeze(sice, ITV)
else
T_fr = 0.0
enth_fp = enth_from_TS(0.0, 0.0, ITV)
endif
max_temp = T_fr ; max_enth = enth_fp
if (present(temp_max)) then ; if (temp_max < T_fr) then
max_temp = temp_max ; max_enth = enth_from_TS(temp_max, sice, ITV)
endif ; endif
enth_in = enth
dtEU = enth_unit * dtt
! Solve m_lay * (enth_new - enth) = dtEU * (htg - fb*t_new)
! t_new = Temp_from_En_S(enth_new, Sice, ITV)
if (m_lay == 0.0) then
new_temp = min(htg / fb, max_temp)
enth = enth_from_TS(new_temp, sice, ITV)
elseif (dtEU * (htg - fb*max_temp) >= m_lay*(max_enth - enth_in)) then
! There is enough heat being applied here that the ice would be above the
! maximum temperature (often the freezing point). The ice should be set to
! the maximum temperature and enthalpy, and the extra heat stored for later
! use in melting.
extra_enth = m_lay*(enth_in - max_enth) + dtEU * (htg - fb*max_temp)
extra_heat = extra_enth / enth_unit
new_temp = max_temp
enth = max_enth
elseif ( sice == 0.0 ) then ! Note that T_fr = 0.
! dT_dEnth is 0 for enth > enth_fp.
! dT_dEnth = dTemp_dEnth_EnS(enth_in, Sice, ITV)
dT_dEnth = 1.0 / (Cp_Ice * enth_unit)
! Solve for enth: m_lay * (enth - enth_in) =
! dtEU * (htg - fb*T_fr - fb*dT_dEnth*(enth - enth_fp))
! enth = enth_in + dtEU * (htg - fb*(0.0 - dT_dEnth*(enth_fp-enth_in))) / &
! (m_lay + dtEU*b*dT_dEnth)
! Or equivalently... (noting that T_fr = 0.0)
enth = enth_fp + (dtEU * (htg - fb*0.0) + m_lay * (enth_in-enth_fp)) / &
(m_lay + dtEU*(fb*dT_dEnth))
! The following is equivalent to new_temp = Temp_from_En_S(enth, 0.0, ITV)
! or new_temp = dT_dEnth * (enth - enth_fp) ! + T_fr==0.
! but it avoids serious roundoff issues later on when b is large.
new_temp = dT_dEnth * ((dtEU * (htg - fb*0.0) + m_lay * (enth_in-enth_fp)) / &
(m_lay + dtEU*(fb*dT_dEnth)))
else! if (Cp_ice == Cp_brine) then
En_J = (enth_in - enthalpy_liquid(0.0, 0.0, ITV)) / enth_unit
! Solve a quadratic equation for the new layer temperature, tn:
!
! m * (En_J - (Cp_Water-Cp_Ice)*T_fr + L + htg*dt/m) =
! (m*Cp_Ice + b*dt) *tn + m*LI*T_fr/tn
!
AA = m_lay *Cp_Ice + fb*dtt
BB = -(m_lay*((En_J - (Cp_Water-Cp_Ice)*T_fr) + LI) + htg*dtt)
CC = m_lay * LI*T_fr
! This form avoids round-off errors.
if (BB >= 0) then
new_temp = -(BB + sqrt(BB*BB - 4*AA*CC)) / (2*AA)
else
new_temp = (2*CC) / (-BB + sqrt(BB*BB - 4*AA*CC))
endif
if (Cp_ice == Cp_brine) then
enth = enth_from_TS(new_temp, sice, ITV)
else ! (Cp_ice /= Cp_brine)
! Correct the new temperature estimate.
! Solve for enth & -273.15 < T_g < T_fr < 0
! m_lay*(enth - En_J) = dtt * (htg - fb*T_g)
! enth = (-LI * (1.0 - T_fr/T_g)) + &
! ((Cp_Ice*Tg + TfxdCP_WI) - TfxdCp_BI*log(T_fr/T_g))
! Err = m_lay*(enth - En_J) + dtt * (fb*T_g - htg)
! Err = m_lay*((-LI * (1.0 - T_fr/T_g)) + &
! ((Cp_Ice*Tg + TfxdCP_WI) - TfxdCp_BI*log(T_fr/T_g)) - En_J) + dt * (fb*T_g - htg)
! En_J = enth_in / enth_unit - enth_liq_0
TfxdCp_WI = T_fr*(Cp_Water-Cp_Ice)
TfxdCp_BI = T_fr*(Cp_Brine-Cp_Ice)
Err_Tind = (m_lay*(-LI + TfxdCP_WI - En_J) - dtt*htg)
T_min = -273.15
Err_Tmin = m_lay*(LI * (T_fr/T_min) + (Cp_Ice*T_min - TfxdCp_BI*log(T_fr/T_min))) + &
(dtt * fb * T_min + Err_Tind)
! Approximate this as
! Err_Tmin = m_lay*((Cp_Ice*T_min - TfxdCp_BI*log(T_fr/T_min))) + &
! (dtt * fb * T_min + Err_Tind) ?
T_max = T_fr
Err_Tmax = m_lay*(T_fr*Cp_Water - En_J) + dtt * (fb * T_fr - htg)
T_g = new_temp
! Using a false position method first-guess instead adds about 2 iterations.
! T_g = (Err_Tmax * T_min - Err_Tmin * T_max) / (Err_Tmax - Err_Tmin)
! T_itt(:) = 0.0 ; dTemp(:) = 0.0 ; Err_itt(:) = 0.0
do itt=1,20 ! Note that 3 or 4 iterations usually are enough.
Err = m_lay*(LI * (T_fr/T_g) + (Cp_Ice*T_g - TfxdCp_BI*log(T_fr/T_g))) + &
(dtt * fb * T_g + Err_Tind)
! T_itt(itt) = T_g ; Err_itt(itt) = Err
if (abs(Err) <= Enth_tol*(abs(Err_Tind) + m_lay*LI + abs(dtt*fb*T_fr))) then
new_temp = T_g ; exit
elseif (Err < 0.0) then
T_min = T_g ; Err_Tmin = Err
else
T_max = T_g ; Err_Tmax = Err
endif
! Use the more efficient Newton's method of McDougall & Witherspoon (2014),
! Appl. Math. Lett., 29, 20-25.
T_deriv = T_g
if (itt > 1) then ! Reuse the estimate of dT_dEn from the last iteration.
if ((dErr_dT*T_g - 0.5*Err > dErr_dT*T_min) .and. &
(dErr_dT*T_g - 0.5*Err < dErr_dT*T_max)) &
T_deriv = T_g - 0.5*Err / dErr_dT
endif
dErr_dT = m_lay*( -LI * (T_fr / T_deriv**2) + &
(Cp_Ice + TfxdCp_BI / T_deriv)) + dtt*fb ! >= 0.0
! T_prev = T_g
if ((dErr_dT*T_g - Err > dErr_dT*T_min) .and. &
(dErr_dT*T_g - Err < dErr_dT*T_max)) then
T_g = T_g - Err / dErr_dT
else
T_g = (Err_Tmax * T_min - Err_Tmin * T_max) / &
(Err_Tmax - Err_Tmin)
endif
! dTemp(itt) = T_g - T_prev
enddo
new_temp = T_g ! Use the best guess.
enth = enth_from_TS(new_temp, sice, ITV)
! write (*,'("T_itt = ",8F14.8)') T_itt(1:8)
! write (*,'("Err_itt = ",8(1PE14.6))') Err_itt(1:8)
! write (*,'("dTemp = ",8(1Pe12.4))') dTemp(1:8)
! if (m_lay > 1e-5) then
! ! Check the answers...
! write (*,'(" Enth, Enth_in, Enth_err = ",3(1Pe14.4))') enth, En_J, &
! (enth - En_J) - dtt * (htg - fb*new_temp) / m_lay
! endif
endif
endif
! Figure out whether it is more accurate to use the explicit form for the
! fluxes or to invert enthalpy conservation for the fluxes.
if (abs(hf_err_rat*dftop_dT) <= m_lay) then
! This branch currently applies in all interior ice layers.
ftop = ftop_in + dftop_dT*new_temp
! An explicit estimate (or no update) is used for the top flux.
if (hf_err_rat*dfbot_dT <= m_lay) then ! Use the explicit expression for fbot.
fbot = fbot_in + dfbot_dT*new_temp
else ! Use conservation to invert for fbot.