forked from NOAA-GFDL/SIS2
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathSIS2_ice_thm.F90
2056 lines (1763 loc) · 92.9 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) !
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
module SIS2_ice_thm
! for calling delta-Eddington shortwave from ice_optics
use ice_shortwave_dEdd, only : shortwave_dEdd0_set_snow, shortwave_dEdd0_set_pond
use ice_shortwave_dEdd, only : shortwave_dEdd0, shortwave_dEdd0_set_params
use ice_shortwave_dEdd, only : dbl_kind, int_kind, nilyr, nslyr
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 constants_mod, only : hlv, hlf ! latent heats of vaporization and fusion.
implicit none ; private
public :: get_thermo_coefs, get_SIS2_thermo_coefs, SIS2_ice_thm_end
public :: SIS2_ice_thm_init, ice_optics_SIS2, ice_temp_SIS2
public :: ice_resize_SIS2, add_frazil_SIS2, rebalance_ice_layers
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_SeaWater ! The heat capacity of liquid seawater, in J/(kg K).
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_SeaWater, but for
! algorithmic convenience has often been
! 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 :: dTf_dS ! The derivative of the freezing point with salinity,
! in degC per PSU. (dTf_dS is negative.)
real :: mu_TS ! Negative the derivative of the freezing point with
! salinity, in degC per PSU.
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.
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)
! albedos are from CSIM4 assumming 0.53 visible and 0.47 near-ir insolation
real :: alb_snow ! albedo of snow (not melting)
real :: alb_ice ! albedo of ice (not melting)
real :: pen_ice ! ice surface penetrating solar fraction
real :: opt_dep_ice ! ice optical depth (m)
real :: t_range_melt ! melt albedos scaled in below melting T
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 :: old_heat_cap ! If true, use an older linearization of the
! sea ice heat capacity for temperatures above the
! freezing pointin the calculation of the initial ice
! temperature estimate. The default is false.
logical :: do_deltaEdd = .true. ! If true, use a delta-Eddington radiative
! transfer calculation for the shortwave radiation
! within the sea-ice and snow.
end type SIS2_ice_thm_CS
contains
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
! ice_thm_param - set ice thermodynamic parameters !
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
subroutine SIS2_ice_thm_init(param_file, CS, ITV, init_EOS )
use constants_mod, only : hlf ! latent heat of fusion - 334e3 J/(kg-ice)
type(param_file_type), intent(in) :: param_file
type(SIS2_ice_thm_CS), pointer :: CS
type(ice_thermo_type), pointer :: ITV ! A pointer to the ice thermodynamic parameter structure.
logical, optional, intent(in) :: init_EOS
real :: sal_ice_top(1) ! A specified surface salinity of ice.
real :: deltaEdd_R_ice ! Mysterious delta-Eddington tuning parameters, unknown.
real :: deltaEdd_R_snow ! Mysterious delta-Eddington tuning parameters, unknown.
real :: deltaEdd_R_pond ! Mysterious delta-Eddington tuning parameters, unknown.
character(len=40) :: mod = "SIS2_ice_thm" ! This module's name.
if (.not.associated(CS)) allocate(CS)
if (.not.associated(ITV)) allocate(ITV)
! LI must be taken from the constants mod for internal consistency in the
! coupled climate model.
ITV%LI = hlf
call get_param(param_file, mod, "RHO_OCEAN", ITV%Rho_water, &
"The nominal density of sea water as used by SIS.", &
units="kg m-3", default=1030.0)
call get_param(param_file, mod, "RHO_ICE", ITV%Rho_ice, &
"The nominal density of sea ice as used by SIS.", &
units="kg m-3", default=905.0)
call get_param(param_file, mod, "RHO_SNOW", ITV%Rho_snow, &
"The nominal density of snow as used by SIS.", &
units="kg m-3", default=330.0)
call get_param(param_file, mod, "CP_ICE", ITV%Cp_ice, &
"The heat capacity of fresh ice, approximated as a \n"//&
"constant.", units="J kg-1 K-1", default=2100.0)
call get_param(param_file, mod, "CP_SEAWATER", ITV%Cp_SeaWater, &
"The heat capacity of sea water, approximated as a \n"//&
"constant.", units="J kg-1 K-1", default=4200.0)
call get_param(param_file, mod, "CP_WATER", ITV%Cp_water, &
"The heat capacity of water in sea-ice, approximated as \n"//&
"a constant. CP_WATER and CP_SEAWATER should be equal, \n"//&
"but for computational convenience CP_WATER has often \n"//&
"been set equal to CP_ICE instead.", units="J kg-1 K-1", &
default=ITV%Cp_SeaWater)
call get_param(param_file, mod, "CP_BRINE", ITV%Cp_brine, &
"The heat capacity of water in brine pockets within the \n"//&
"sea-ice, approximated as a constant. CP_BRINE and \n"//&
"CP_WATER should be equal, but for computational \n"//&
"convenience CP_BRINE has often been set equal to CP_ICE.", &
units="J kg-1 K-1", default=ITV%Cp_ice) !### CHANGE LATER TO default=CP_WATER)
call get_param(param_file, mod, "DTFREEZE_DS", ITV%dTf_dS, &
"The derivative of the freezing temperature with salinity.", &
units="deg C PSU-1", default=-0.054)
ITV%mu_TS = -ITV%dTf_dS
call get_thermo_coefs(ice_salinity=sal_ice_top)
CS%temp_ice_freeze = T_freeze(sal_ice_top(1), ITV)
call get_param(param_file, mod, "ENTHALPY_LIQUID_0", ITV%enth_liq_0, &
"The enthalpy of liquid fresh water at 0 C. The solutions \n"//&
"should be physically consistent when this is adjusted, \n"//&
"because only the relative value is of physical meaning, \n"//&
"but roundoff errors can change the solution.", units="J kg-1", &
default=0.0)
call get_param(param_file, mod, "ENTHALPY_UNITS", ITV%enth_unit, &
"A constant that rescales enthalpy from J/kg to a \n"//&
"different scale in its internal representation. Changing \n"//&
"this by a power of 2 is useful for debugging, as answers \n"//&
"should not change. A negative values is taken as an inverse.", &
units="J kg-1", default=1.0)
if (ITV%enth_unit < 0.) ITV%enth_unit = -1.0 / ITV%enth_unit
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_DELTA_EDDINGTON_SW", CS%do_deltaEdd, &
"If true, a delta-Eddington radiative transfer calculation \n"//&
"for the shortwave radiation within the sea-ice.", default=.true.)
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 get_param(param_file, mod, "OLD_ICE_HEAT_CAPACITY", CS%old_heat_cap, &
"If true, use an older linearization of the sea ice heat \n"//&
"capacity for temperatures above the freezing point in the \n"//&
"calculation of the initial ice temperature estimate.", default=.false.)
if (CS%do_deltaEdd) then
call get_param(param_file, mod, "ICE_DELTA_EDD_R_ICE", deltaEdd_R_ice, &
"A dreadfully documented tuning parameter for the radiative \n"//&
"propeties of sea ice with the delta-Eddington radiative \n"//&
"transfer calculation.", units="perhaps nondimensional?", default=0.0)
call get_param(param_file, mod, "ICE_DELTA_EDD_R_SNOW", deltaEdd_R_snow, &
"A dreadfully documented tuning parameter for the radiative \n"//&
"propeties of snow on sea ice with the delta-Eddington \n"//&
"radiative transfer calculation.", &
units="perhaps nondimensional?", default=0.0)
call get_param(param_file, mod, "ICE_DELTA_EDD_R_POND", deltaEdd_R_pond, &
"A dreadfully documented tuning parameter for the radiative \n"//&
"propeties of meltwater ponds on sea ice with the delta-Eddington \n"//&
"radiative transfer calculation.", units="perhaps nondimensional?", &
default=0.0)
call shortwave_dEdd0_set_params(deltaEdd_R_ice, deltaEdd_R_snow, deltaEdd_R_pond)
else
call get_param(param_file, mod, "SNOW_ALBEDO", CS%alb_snow, &
"The albedo of dry snow atop sea ice.", units="nondim", &
default=0.85)
call get_param(param_file, mod, "ICE_ALBEDO", CS%alb_ice, &
"The albedo of dry bare sea ice.", units="nondim", &
default=0.5826)
call get_param(param_file, mod, "ICE_SW_PEN_FRAC", CS%pen_ice, &
"The fraction of the unreflected shortwave radiation that \n"//&
"penetrates into the ice.", units="Nondimensional", default=0.3)
call get_param(param_file, mod, "ICE_OPTICAL_DEPTH", CS%opt_dep_ice, &
"The optical depth of shortwave radiation in sea ice.", &
units="m", default=0.67)
call get_param(param_file, mod, "ALBEDO_T_MELT_RANGE", CS%t_range_melt, &
"The temperature range below freezing over which the \n"//&
"albedos are changed by partial melting.", units="degC", &
default=1.0)
endif
if (present(init_EOS)) then ; if (init_EOS) then
if (.not.associated(ITV%EOS)) call EOS_init(param_file, ITV%EOS)
endif ; endif
end subroutine SIS2_ice_thm_init
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
! ice_optics - set albedo, penetrating solar, and ice/snow transmissivity !
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
subroutine ice_optics_SIS2(hs, hi, ts, tfw, NkIce, alb_vis_dir, alb_vis_dif, &
alb_nir_dir, alb_nir_dif, abs_sfc, abs_snow, abs_ice_lay, &
abs_ocn, abs_int, CS, coszen_in)
real, intent(in ) :: hs ! snow thickness (m-snow)
real, intent(in ) :: hi ! ice thickness (m-ice)
real, intent(in ) :: ts ! surface temperature
real, intent(in ) :: tfw ! seawater freezing temperature
integer, intent(in) :: NkIce
real, intent( out) :: alb_vis_dir ! ice surface albedo (0-1)
real, intent( out) :: alb_vis_dif ! ice surface albedo (0-1)
real, intent( out) :: alb_nir_dir ! ice surface albedo (0-1)
real, intent( out) :: alb_nir_dif ! ice surface albedo (0-1)
real, intent( out) :: abs_sfc ! frac abs sw abs at surface
real, intent( out) :: abs_snow ! frac abs sw abs in snow
real, intent( out) :: abs_ice_lay(NkIce) ! frac abs sw abs by each ice layer
real, intent( out) :: abs_ocn ! frac abs sw abs in ocean
real, intent( out) :: abs_int ! frac abs sw abs in ice interior
type(SIS2_ice_thm_CS), intent(in) :: CS
real, intent(in),optional :: coszen_in
real :: alb, as, ai, snow_cover, fh
real :: coalb, I_coalb ! The coalbedo and its reciprocal.
real :: SW_frac_top ! The fraction of the SW at the top of the snow that
! is still present at the top of each ice layer (ND).
real :: opt_decay_lay ! The optical extinction in each ice layer (ND).
real :: pen ! frac sw passed below the surface (frac 1-pen absorbed at the surface)
integer :: m
character(len=200) :: mesg
integer (kind=int_kind) :: &
nx_block, ny_block, & ! block dimensions
icells ! number of ice-covered grid cells
integer (kind=int_kind), dimension (1) :: &
indxi , & ! compressed indices for ice-covered cells
indxj
! inputs
real (kind=dbl_kind), dimension (1,1) :: &
aice , & ! concentration of ice
vice , & ! volume of ice
vsno , & ! volume of snow
Tsfc , & ! surface temperature
coszen , & ! cosine of solar zenith angle
tarea , & ! cell area - not used
swvdr , & ! sw down, visible, direct (W/m^2)
swvdf , & ! sw down, visible, diffuse (W/m^2)
swidr , & ! sw down, near IR, direct (W/m^2)
swidf ! sw down, near IR, diffuse (W/m^2)
! outputs
real (kind=dbl_kind), dimension (1,1) :: &
fs , & ! horizontal coverage of snow
fp , & ! pond fractional coverage (0 to 1)
hp ! pond depth (m)
real (kind=dbl_kind), dimension (1,1,1) :: &
rhosnw , & ! density in snow layer (kg/m3)
rsnw ! grain radius in snow layer (micro-meters)
real (kind=dbl_kind), dimension (1,1,18) :: &
trcr ! aerosol tracers
real (kind=dbl_kind), dimension (1,1) :: &
alvdr , & ! visible, direct, albedo (fraction)
alvdf , & ! visible, diffuse, albedo (fraction)
alidr , & ! near-ir, direct, albedo (fraction)
alidf , & ! near-ir, diffuse, albedo (fraction)
fswsfc , & ! SW absorbed at snow/bare ice/pondedi ice surface (W m-2)
fswint , & ! SW interior absorption (below surface, above ocean,W m-2)
fswthru ! SW through snow/bare ice/ponded ice into ocean (W m-2)
real (kind=dbl_kind), dimension (1,1,1) :: &
Sswabs ! SW absorbed in snow layer (W m-2)
real (kind=dbl_kind), dimension (1,1,NkIce) :: &
Iswabs ! SW absorbed in ice layer (W m-2)
real (kind=dbl_kind), dimension (1,1) :: &
albice , & ! bare ice albedo, for history
albsno , & ! snow albedo, for history
albpnd ! pond albedo, for history
if (CS%do_deltaEdd) then
if (nilyr /= NkIce) then
write(mesg, '("The Delta-Eddington sea-ice radiation is hard-coded to use ",(I4),&
&" ice layers, not ",(I4),".")') nilyr, NkIce
call SIS_error(FATAL, mesg)
endif
! temporary for delta-Eddington shortwave call
nx_block = 1 ; ny_block = 1
icells = 1 ; indxi(1) = 1 ; indxj(1) = 1
aice(1,1) = 1.0
tarea(1,1) = 1.0 ! not used
! stuff that matters
coszen(1,1) = cos(3.14*67.0/180.0) ! NP summer solstice
if(present(coszen_in)) coszen(1,1) = max(0.01,coszen_in)
Tsfc(1,1) = ts
vsno(1,1) = hs
vice(1,1) = hi
swvdr(1,1) = 0.25
swvdf(1,1) = 0.25
swidr(1,1) = 0.25
swidf(1,1) = 0.25
call shortwave_dEdd0_set_snow(nx_block, ny_block, icells, indxi, indxj, &
aice, vsno, Tsfc, fs, rhosnw, rsnw) ! out: fs, rhosnw, rsnw
call shortwave_dEdd0_set_pond(nx_block, ny_block, icells, indxi, indxj, &
aice, Tsfc, fs, fp, hp) ! out: fp, hp
call shortwave_dEdd0 (nx_block, ny_block, icells, indxi, indxj, coszen, &
aice, vice, vsno, fs, rhosnw, rsnw, fp, hp, swvdr, swvdf, &
swidr, swidf, alvdf, alvdr, alidr, alidf, fswsfc, fswint, &
fswthru, Sswabs, Iswabs, albice, albsno, albpnd)
! out: alvdf, alvdr, and subsequent.
! Note: fswint = Sswabs + sum(Iswabs)
alb = 1.0 - (fswsfc(1,1) + (fswint(1,1) + fswthru(1,1)))
coalb = fswsfc(1,1) + (fswint(1,1) + fswthru(1,1))
I_coalb = 0.0 ; if (coalb > 0.0) I_coalb = 1.0 / coalb
abs_sfc = fswsfc(1,1) * I_coalb
abs_snow = Sswabs(1,1,1) * I_coalb
do m=1,NkIce ; abs_ice_lay(m) = Iswabs(1,1,m) * I_coalb ; enddo
abs_ocn = fswthru(1,1) * I_coalb
alb_vis_dir = alvdr(1,1)
alb_vis_dif = alvdf(1,1)
alb_nir_dir = alidr(1,1)
alb_nir_dif = alidf(1,1)
pen = (fswint(1,1) + fswthru(1,1)) * I_coalb
abs_int = fswint(1,1) * I_coalb
else
as = CS%alb_snow ; ai = CS%alb_ice
snow_cover = hs/(hs+0.02) ! thin snow partially covers ice
fh = min(atan(5.0*hi)/atan(5.0*0.5),1.0) ! use this form from CSIM4 to
! reduce albedo for thin ice
if (ts+CS%T_RANGE_MELT > CS%temp_ice_freeze) then ! reduce albedo for melting as in
! CSIM4 assuming 0.53/0.47 vis/ir
as = as-0.1235*min((ts+CS%T_RANGE_MELT-CS%temp_ice_freeze)/CS%T_RANGE_MELT,1.0)
ai = ai-0.075 *min((ts+CS%T_RANGE_MELT-CS%temp_ice_freeze)/CS%T_RANGE_MELT,1.0)
endif
ai = fh*ai+(1-fh)*0.06 ! reduce albedo for thin ice
alb = snow_cover*as + (1-snow_cover)*ai
alb_vis_dir = alb ; alb_vis_dif = alb
alb_nir_dir = alb ; alb_nir_dif = alb
pen = (1-snow_cover)*CS%pen_ice
opt_decay_lay = exp(-hi/(NkIce*CS%opt_dep_ice))
abs_ocn = pen * exp(-hi/CS%opt_dep_ice)
abs_sfc = 1.0 - pen
abs_snow = 0.0
SW_frac_top = pen
do m=1,NkIce
abs_ice_lay(m) = SW_frac_top * (1.0 - opt_decay_lay)
SW_frac_top = SW_frac_top * opt_decay_lay
enddo
abs_int = pen - SW_frac_top
!! check for ice albedos out of range (0 to 1)
! if (alb.lt.0.0 .or. alb.gt.1.0) then
! print *,'ice_optics: albedo out of range, alb=',alb
! print *,'cs=',cs, 'as=',as, 'ai=',ai
! print *,'ts=',ts, 'fh=',fh, 'hs=',hs, 'hi=',hi, 'tfw=',tfw
! print *,'ALB_SNO=',ALB_SNO, 'ALB_ICE=',ALB_ICE, 'T_RANGE_MELT,=',T_RANGE_MELT, 'TFI=',TFI
! stop
! end if
endif
end subroutine ice_optics_SIS2
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
! ice_temp_SIS2 - A subroutine that calculates the snow and ice enthalpy !
! changes due to surface forcing. !
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
subroutine ice_temp_SIS2(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_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 :: I_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)
A = -sh_T0
I_enth_unit = 1.0 / ITV%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 / ITV%Rho_ice, CS%H_LO_LIM)
hsnow_eff = mL_snow / ITV%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*ITV%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 * ITV%Cp_ice
elseif ((temp_IC(k) < tfi(k)) .or. (CS%old_heat_cap)) then
bb(k) = mL_ice * (ITV%Cp_ice - (tfi(k) / temp_IC(k)**2) * &
(ITV%LI - (ITV%Cp_brine-ITV%Cp_ice) * temp_IC(k)) )
! Or mroe generally:
! S_Sf = sice(k) / calculate_S_freeze(temp_IC(k))
! bb(k) = mL_ice * (ITV%Cp_ice + dSf_dT*ITV%LI + S_Sf * &
! ((ITV%Cp_brine-ITV%Cp_ice))
else ! Use the value when temp_IC = tfi.
bb(k) = mL_ice * (ITV%Cp_brine - ITV%LI / 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) then
! 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 * ITV%Cp_Ice + ITV%LI)
e_extra = 0.0
if (tsurf > tsf) then ! The surface is melting: update enthalpy with the surface at melt temp.
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*ITV%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*ITV%enth_unit)
if (abs(d_e_extra) > 1.0e-12*(abs(col_enth3) + abs(col_enth2b) + &
abs(e_extra_sum*ITV%enth_unit))) then
d_e_extra = (col_enth3 - col_enth2b) - e_extra_sum*ITV%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, LI
integer :: itt
! real :: T_itt(20), dTemp(20), Err_itt(20)
Cp_Ice = ITV%Cp_Ice ; LI = ITV%LI
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 = ITV%Cp_Water*(tp - T_fr) ! >= 0
else
E0 = Cp_Ice*(tp - T_fr) - LI*(1 - T_fr/tp) ! < 0
if (ITV%Cp_ice /= ITV%Cp_brine) E0 = E0 + (ITV%Cp_brine - ITV%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)) / (ITV%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 (ITV%Cp_ice /= ITV%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*(ITV%Cp_Brine-ITV%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.