forked from NOAA-GFDL/SIS2
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathSIS_slow_thermo.F90
1338 lines (1174 loc) · 65.9 KB
/
SIS_slow_thermo.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
!***********************************************************************
!* GNU General Public License *
!* This file is a part of SIS2. *
!* *
!* SIS2 is free software; you can redistribute it and/or modify it and *
!* are expected to follow the terms of the GNU General Public License *
!* as published by the Free Software Foundation; either version 2 of *
!* the License, or (at your option) any later version. *
!* *
!* SIS2 is distributed in the hope that it will be useful, but WITHOUT *
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY *
!* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public *
!* License for more details. *
!* *
!* For the full text of the GNU General Public License, *
!* write to: Free Software Foundation, Inc., *
!* 675 Mass Ave, Cambridge, MA 02139, USA. *
!* or see: http://www.gnu.org/licenses/gpl.html *
!***********************************************************************
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
! SIS2 is a SEA ICE MODEL for coupling through the GFDL exchange grid. SIS2 !
! is a revision of the original SIS with have extended capabilities, including !
! the option of using a B-grid or C-grid spatial discretization. The SIS2 !
! software has been extensively reformulated from SIS for greater consistency !
! with the Modular Ocean Model, version 6 (MOM6), and to permit might tighter !
! dynamical coupling between the ocean and sea-ice. !
! This module handles the main updates of the ice states at the slower time- !
! scales of the couplng or the interactions with the ocean, including the ice !
! mass balance and related thermodynamics and salinity changes, and !
! thermodynamic coupling with the ocean. The radiative heating and diffusive !
! temperature changes due to coupling with the atmosphere are handled elsewhere. !
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
module SIS_slow_thermo
use SIS_diag_mediator, only : enable_SIS_averaging, disable_SIS_averaging
use SIS_diag_mediator, only : post_SIS_data, post_data=>post_SIS_data
use SIS_diag_mediator, only : query_SIS_averaging_enabled, SIS_diag_ctrl
use SIS_diag_mediator, only : register_diag_field=>register_SIS_diag_field
use SIS_sum_output, only : SIS_sum_out_CS, write_ice_statistics! , SIS_sum_output_init
use SIS_sum_output, only : accumulate_bottom_input, accumulate_input_1, accumulate_input_2
! use MOM_domains, only : pass_var
! ! use MOM_dyn_horgrid, only : dyn_horgrid_type, create_dyn_horgrid, destroy_dyn_horgrid
use MOM_error_handler, only : SIS_error=>MOM_error, FATAL, WARNING, SIS_mesg=>MOM_mesg
use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint
use MOM_file_parser, only : get_param, log_param, log_version, param_file_type
use MOM_hor_index, only : hor_index_type
use MOM_EOS, only : EOS_type, calculate_density_derivs
use fms_mod, only : clock_flag_default
use mpp_mod, only : mpp_clock_id, mpp_clock_begin, mpp_clock_end
use mpp_mod, only : CLOCK_COMPONENT, CLOCK_LOOP, CLOCK_ROUTINE
use time_manager_mod, only : time_type, time_type_to_real! , get_date, get_time
! use time_manager_mod, only : set_date, set_time, operator(+), operator(-)
! use MOM_time_manager, only : operator(>), operator(*), operator(/), operator(/=)
use data_override_mod, only : data_override
use SIS_types, only : ice_state_type, ice_ocean_flux_type, fast_ice_avg_type
use SIS_types, only : ocean_sfc_state_type
use SIS_types, only : IST_chksum, IST_bounds_check
use ice_utils_mod, only : post_avg
use SIS_hor_grid, only : SIS_hor_grid_type
use ice_grid, only : ice_grid_type
use ice_spec_mod, only : get_sea_surface
use SIS_tracer_flow_control, only : SIS_call_tracer_column_fns
use SIS2_ice_thm, only : SIS2_ice_thm_CS, SIS2_ice_thm_init, SIS2_ice_thm_end
use SIS2_ice_thm, only : ice_resize_SIS2, add_frazil_SIS2, rebalance_ice_layers
use SIS2_ice_thm, only : get_SIS2_thermo_coefs, enthalpy_liquid_freeze
use SIS2_ice_thm, only : enth_from_TS, Temp_from_En_S
use SIS2_ice_thm, only : T_freeze, enthalpy_liquid, calculate_T_freeze
use ice_transport_mod, only : adjust_ice_categories, ice_transport_CS
use SIS_tracer_flow_control, only : SIS_tracer_flow_control_CS
implicit none ; private
#include <SIS2_memory.h>
public :: slow_thermodynamics, SIS_slow_thermo_init, SIS_slow_thermo_end
public :: slow_thermo_CS, SIS_slow_thermo_set_ptrs
type slow_thermo_CS ; private
logical :: specified_ice ! If true, the sea ice is specified and there is
! no need for ice dynamics.
real :: ice_bulk_salin ! The globally constant sea ice bulk salinity, in g/kg
! that is used to calculate the ocean salt flux.
real :: ice_rel_salin ! The initial bulk salinity of sea-ice relative to the
! salinity of the water from which it formed, nondim.
logical :: filling_frazil ! If true, apply frazil to fill as many categories
! as possible to fill in a uniform (minimum) amount
! of frazil in all the thinnest categories.
! Otherwise the frazil is always assigned to a
! single category with part size > 0.01.
real :: fraz_fill_time ! A timescale with which the filling frazil causes
! the thinest cells to attain similar thicknesses,
! or a negative number to apply the frazil flux
! uniformly, in s.
logical :: do_ridging ! If true, apply a ridging scheme to the convergent
! ice. The original SIS2 implementation is based on
! work by Torge Martin. Otherwise, ice is compressed
! proportionately if the concentration exceeds 1.
logical :: do_ice_restore ! If true, restore the sea-ice toward climatology
! by applying a restorative heat flux.
real :: ice_restore_timescale ! The time scale for restoring ice when
! do_ice_restore is true, in days.
logical :: do_ice_limit ! Limit the sea ice thickness to max_ice_limit.
real :: max_ice_limit ! The maximum sea ice thickness, in m, when
! do_ice_limit is true.
logical :: nudge_sea_ice = .false. ! If true, nudge sea ice concentrations towards observations.
real :: nudge_sea_ice_rate = 0.0 ! The rate of cooling of ice-free water that
! should be ice covered in order to constrained the
! ice concentration to track observations. A suggested
! value is of order 10000 W m-2.
real :: nudge_stab_fac ! A factor that determines whether the buoyancy
! flux associated with the sea ice nudging of
! warm water includes a freshwater flux so as to
! be destabilizing on net (<1), stabilizing (>1),
! or neutral (=1). The default is 1.
real :: nudge_conc_tol ! The tolerance for mismatch in the sea ice concentations
! before nudging begins to be applied.
logical :: debug ! If true, write verbose checksums for debugging purposes.
logical :: column_check ! If true, enable the heat check column by column.
real :: imb_tol ! The tolerance for imbalances to be flagged by
! column_check, nondim.
logical :: bounds_check ! If true, check for sensible values of thicknesses
! temperatures, fluxes, etc.
integer :: n_calls = 0 ! The number of times update_ice_model_slow_down
! has been called.
type(time_type), pointer :: Time ! A pointer to the ocean model's clock.
type(SIS_diag_ctrl), pointer :: diag ! A structure that is used to regulate the
! timing of diagnostic output.
! These are pointers to the control structures for subsidiary modules.
type(SIS2_ice_thm_CS), pointer :: ice_thm_CSp => NULL()
type(ice_transport_CS), pointer :: ice_transport_CSp => NULL()
type(SIS_sum_out_CS), pointer :: sum_output_CSp => NULL()
type(SIS_tracer_flow_control_CS), pointer :: tracer_flow_CSp => NULL()
integer :: id_qflim=-1, id_qfres=-1, id_fwnudge=-1
integer :: id_lsrc=-1, id_lsnk=-1, id_bsnk=-1, id_sn2ic=-1
end type slow_thermo_CS
integer :: iceClock5, iceClock6, iceClock7
contains
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
! post_flux_diagnostics - write out any diagnostics of surface fluxes. !
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
subroutine post_flux_diagnostics(IST, FIA, IOF, CS, G, IG, Idt_slow)
type(ice_state_type), intent(in) :: IST
type(fast_ice_avg_type), intent(in) :: FIA
type(ice_ocean_flux_type), intent(in) :: IOF
type(slow_thermo_CS), pointer :: CS
type(SIS_hor_grid_type), intent(in) :: G
type(ice_grid_type), intent(in) :: IG
real, intent(in) :: Idt_slow
real, dimension(G%isd:G%ied,G%jsd:G%jed) :: tmp2d, net_sw
integer :: i, j, k, m, n, isc, iec, jsc, jec, ncat
isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; ncat = IG%CatIce
! Flux diagnostics
!
if (FIA%id_runoff>0) &
call post_data(FIA%id_runoff, FIA%runoff, CS%diag)
if (FIA%id_calving>0) &
call post_data(FIA%id_calving, FIA%calving_preberg, CS%diag)
if (FIA%id_runoff_hflx>0) &
call post_data(FIA%id_runoff_hflx, FIA%runoff_hflx, CS%diag)
if (FIA%id_calving_hflx>0) &
call post_data(FIA%id_calving_hflx, FIA%calving_hflx_preberg, CS%diag)
! The frazil diagnostic is with the other ocean surface diagnostics.
! if (IST%id_frazil>0) &
! call post_data(IST%id_frazil, FIA%frazil_left*Idt_slow, CS%diag)
if (FIA%id_sh>0) call post_avg(FIA%id_sh, FIA%flux_t_top, IST%part_size, &
CS%diag, G=G)
if (FIA%id_lh>0) call post_avg(FIA%id_lh, FIA%flux_lh_top, IST%part_size, &
CS%diag, G=G)
if (FIA%id_evap>0) call post_avg(FIA%id_evap, FIA%flux_q_top, IST%part_size, &
CS%diag, G=G)
if (FIA%id_slp>0) &
call post_data(FIA%id_slp, FIA%p_atm_surf, CS%diag)
if ((FIA%id_sw>0) .or. (FIA%id_albedo>0)) then
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,ncat,net_sw,IST,FIA)
do j=jsc,jec
do i=isc,iec ; net_sw(i,j) = 0.0 ; enddo
do k=0,ncat ; do i=isc,iec
net_sw(i,j) = net_sw(i,j) + IST%part_size(i,j,k) * ( &
FIA%flux_sw_vis_dir_top(i,j,k) + FIA%flux_sw_vis_dif_top(i,j,k) + &
FIA%flux_sw_nir_dir_top(i,j,k) + FIA%flux_sw_nir_dif_top(i,j,k) )
enddo ; enddo
enddo
if (FIA%id_sw>0) call post_data(FIA%id_sw, net_sw, CS%diag)
if (FIA%id_albedo>0) then
do j=jsc,jec ; do i=isc,iec
if (G%mask2dT(i,j)<=0.5) then
tmp2d(i,j) = -1.0 ! This is land.
elseif ((FIA%flux_sw_dn(i,j) > 0.0)) then
! The 10.0 below is deliberate. An albedo of down to -9 can be reported
! for detecting inconsistent net_sw and sw_dn.
tmp2d(i,j) = (FIA%flux_sw_dn(i,j) - min(net_sw(i,j), 10.0*FIA%flux_sw_dn(i,j))) / &
FIA%flux_sw_dn(i,j)
else
tmp2d(i,j) = 0.0 ! What does the albedo mean at night?
endif
enddo ; enddo
call post_data(FIA%id_albedo, tmp2d, CS%diag)
endif
endif
if (FIA%id_lw>0) call post_avg(FIA%id_lw, FIA%flux_lw_top, &
IST%part_size, CS%diag, G=G)
if (FIA%id_snofl>0) call post_avg(FIA%id_snofl, FIA%fprec_top, &
IST%part_size, CS%diag, G=G)
if (FIA%id_rain>0) call post_avg(FIA%id_rain, FIA%lprec_top, &
IST%part_size, CS%diag, G=G)
if (FIA%id_sw_dn>0) call post_data(FIA%id_sw_dn, FIA%flux_sw_dn, CS%diag)
if (FIA%id_sw_vis>0) then
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,ncat,tmp2d,IST,FIA)
do j=jsc,jec
do i=isc,iec ; tmp2d(i,j) = 0.0 ; enddo
do k=0,ncat ; do i=isc,iec
tmp2d(i,j) = tmp2d(i,j) + IST%part_size(i,j,k) * ( &
FIA%flux_sw_vis_dir_top(i,j,k) + FIA%flux_sw_vis_dif_top(i,j,k) )
enddo ; enddo
enddo
call post_data(FIA%id_sw_vis, tmp2d, CS%diag)
endif
if (FIA%id_sw_dir>0) then
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,ncat,tmp2d,IST,FIA)
do j=jsc,jec
do i=isc,iec ; tmp2d(i,j) = 0.0 ; enddo
do k=0,ncat ; do i=isc,iec
tmp2d(i,j) = tmp2d(i,j) + IST%part_size(i,j,k) * ( &
FIA%flux_sw_vis_dir_top(i,j,k) + FIA%flux_sw_nir_dir_top(i,j,k) )
enddo ; enddo
enddo
call post_data(FIA%id_sw_dir, tmp2d, CS%diag)
endif
if (FIA%id_sw_dif>0) then
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,ncat,tmp2d,IST,FIA)
do j=jsc,jec
do i=isc,iec ; tmp2d(i,j) = 0.0 ; enddo
do k=0,ncat ; do i=isc,iec
tmp2d(i,j) = tmp2d(i,j) + IST%part_size(i,j,k) * ( &
FIA%flux_sw_vis_dif_top(i,j,k) + FIA%flux_sw_nir_dif_top(i,j,k) )
enddo ; enddo
enddo
call post_data(FIA%id_sw_dif, tmp2d, CS%diag)
endif
if (FIA%id_sw_nir_dir>0) call post_avg(FIA%id_sw_nir_dir, FIA%flux_sw_nir_dir_top, &
IST%part_size, CS%diag, G=G)
if (FIA%id_sw_nir_dif>0) call post_avg(FIA%id_sw_nir_dif, FIA%flux_sw_nir_dif_top, &
IST%part_size, CS%diag, G=G)
if (FIA%id_sw_vis_dir>0) call post_avg(FIA%id_sw_vis_dir, FIA%flux_sw_vis_dir_top, &
IST%part_size, CS%diag, G=G)
if (FIA%id_sw_vis_dif>0) call post_avg(FIA%id_sw_vis_dif, FIA%flux_sw_vis_dif_top, &
IST%part_size, CS%diag, G=G)
if (CS%nudge_sea_ice .and. CS%id_fwnudge>0) then
call post_data(CS%id_fwnudge, IOF%melt_nudge, CS%diag)
endif
end subroutine post_flux_diagnostics
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
!> slow_thermodynamics takes care of slow ice thermodynamics and mass changes
subroutine slow_thermodynamics(IST, dt_slow, CS, OSS, FIA, IOF, G, IG)
type(ice_state_type), intent(inout) :: IST
real, intent(in) :: dt_slow ! The thermodynamic step, in s.
type(slow_thermo_CS), pointer :: CS
type(ocean_sfc_state_type), intent(in) :: OSS
type(fast_ice_avg_type), intent(inout) :: FIA
type(ice_ocean_flux_type), intent(inout) :: IOF
type(SIS_hor_grid_type), intent(inout) :: G
type(ice_grid_type), intent(inout) :: IG
real, dimension(SZI_(G),SZJ_(G)) :: &
h_ice_input ! The specified ice thickness, with specified_ice, in m.
real :: rho_ice ! The nominal density of sea ice in kg m-3.
real :: Idt_slow
integer :: i, j, k, l, m, isc, iec, jsc, jec, ncat, NkIce
integer :: isd, ied, jsd, jed
real, dimension(SZI_(G),SZJ_(G),IG%CatIce) :: &
rdg_frac, & ! fraction of ridged ice per category
mi_old ! Ice mass per unit area before thermodynamics.
real :: tmp3 ! This is a bad name - make it more descriptive!
mi_old(:,:,:) = 0.0
isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; ncat = IG%CatIce
isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; NkIce = IG%NkIce
! I_Nk = 1.0 / NkIce
Idt_slow = 0.0 ; if (dt_slow > 0.0) Idt_slow = 1.0/dt_slow
CS%n_calls = CS%n_calls + 1
if (CS%debug) then
call IST_chksum("Start update_ice_model_slow", IST, G, IG)
endif
if (CS%bounds_check) &
call IST_bounds_check(IST, G, IG, "Start of SIS_slow_thermo", OSS=OSS)
! Set the frazil heat flux that remains to be applied. This might need
! to be moved earlier in the algorithm, if there ever to be multiple calls to
! slow_thermodynamics per coupling timestep.
do j=jsc,jec ; do i=isc,iec
FIA%frazil_left(i,j) = OSS%frazil(i,j)
enddo ; enddo
!
! conservation checks: top fluxes
!
call mpp_clock_begin(iceClock7)
call accumulate_input_1(IST, FIA, dt_slow, G, IG, CS%sum_output_CSp)
if (CS%column_check) &
call write_ice_statistics(IST, CS%Time, CS%n_calls, G, IG, CS%sum_output_CSp, &
message=" SIS_slow_thermo", check_column=.true.)
call mpp_clock_end(iceClock7)
!
! Thermodynamics
!
if (CS%specified_ice) then ! over-write changes with specifications.
h_ice_input(:,:) = 0.0
call get_sea_surface(CS%Time, IST%t_surf(isc:iec,jsc:jec,0), IST%part_size(isc:iec,jsc:jec,:), &
h_ice_input(isc:iec,jsc:jec))
call get_SIS2_thermo_coefs(IST%ITV, rho_ice=rho_ice)
do j=jsc,jec ; do i=isc,iec
IST%mH_ice(i,j,1) = h_ice_input(i,j) * (IG%kg_m2_to_H * rho_ice)
enddo ; enddo
endif
! IOF must be updated regardless of whether the ice is specified or the prognostic model
! is being used
if (FIA%num_tr_fluxes>0) then
!It is necessary and sufficient that only one OMP thread goes through the following block
!since IOF is shared between the threads (hence the block is not thread-safe).
!$OMP SINGLE
if (IOF%num_tr_fluxes < 0) then
! This is the first call, and the IOF arrays need to be allocated.
IOF%num_tr_fluxes = FIA%num_tr_fluxes
allocate(IOF%tr_flux_ocn_top(SZI_(G), SZJ_(G), IOF%num_tr_fluxes))
IOF%tr_flux_ocn_top(:,:,:) = 0.0
allocate(IOF%tr_flux_index(size(FIA%tr_flux_index,1), size(FIA%tr_flux_index,2)))
IOF%tr_flux_index(:,:) = FIA%tr_flux_index(:,:)
endif
!$OMP END SINGLE
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,ncat,IST,FIA,IOF)
do m=1,FIA%num_tr_fluxes
do j=jsc,jec ; do i=isc,iec
IOF%tr_flux_ocn_top(i,j,m) = IST%part_size(i,j,0) * FIA%tr_flux_top(i,j,0,m)
enddo ; enddo
do k=1,ncat ; do j=jsc,jec ; do i=isc,iec
IOF%tr_flux_ocn_top(i,j,m) = IOF%tr_flux_ocn_top(i,j,m) + &
IST%part_size(i,j,k) * FIA%tr_flux_top(i,j,k,m)
enddo ; enddo ; enddo
enddo
endif
! No other thermodynamics need to be done for ice that is specified,
if(CS%specified_ice) return ;
! Otherwise, Continue with the remainder of the prognostic slow thermodynamics
!TOM> Store old ice mass per unit area for calculating partial ice growth.
mi_old = IST%mH_ice
!TOM> derive ridged ice fraction prior to thermodynamic changes of ice thickness
! in order to subtract ice melt proportionally from ridged ice volume (see below)
if (CS%do_ridging) then
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,ncat,IST,rdg_frac) &
!$OMP private(tmp3)
do j=jsc,jec ; do k=1,ncat ; do i=isc,iec
tmp3 = IST%mH_ice(i,j,k)*IST%part_size(i,j,k)
rdg_frac(i,j,k) = 0.0 ; if (tmp3 > 0.0) &
rdg_frac(i,j,k) = IST%rdg_mice(i,j,k) / tmp3
enddo ; enddo ; enddo
endif
call enable_SIS_averaging(dt_slow, CS%Time, CS%diag)
! Save out diagnostics of fluxes. This must go before SIS2_thermodynamics.
call post_flux_diagnostics(IST, FIA, IOF, CS, G, IG, Idt_slow)
call disable_SIS_averaging(CS%diag)
call accumulate_input_2(IST, FIA, IOF, IST%part_size, dt_slow, G, IG, CS%sum_output_CSp)
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,IOF)
do j=jsc,jec ; do i=isc,iec
IOF%Enth_Mass_in_atm(i,j) = 0.0 ; IOF%Enth_Mass_out_atm(i,j) = 0.0
IOF%Enth_Mass_in_ocn(i,j) = 0.0 ; IOF%Enth_Mass_out_ocn(i,j) = 0.0
enddo ; enddo
! The thermodynamics routines return updated values of the ice and snow
! masses-per-unit area and enthalpies.
call SIS2_thermodynamics(IST, dt_slow, CS, OSS, FIA, IOF, G, IG)
!TOM> calculate partial ice growth for ridging and aging.
if (CS%do_ridging) then
! ice growth (IST%mH_ice > mi_old) does not affect ridged ice volume
! ice melt (IST%mH_ice < mi_old) reduces ridged ice volume proportionally
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,ncat,IST,mi_old,rdg_frac)
do j=jsc,jec ; do k=1,ncat ; do i=isc,iec
if (IST%mH_ice(i,j,k) < mi_old(i,j,k)) &
IST%rdg_mice(i,j,k) = IST%rdg_mice(i,j,k) + rdg_frac(i,j,k) * &
(IST%mH_ice(i,j,k) - mi_old(i,j,k)) * IST%part_size(i,j,k)
IST%rdg_mice(i,j,k) = max(IST%rdg_mice(i,j,k), 0.0)
enddo ; enddo ; enddo
endif
! Other routines that do thermodynamic vertical processes should be added here
! Do tracer column physics
call enable_SIS_averaging(dt_slow, CS%Time, CS%diag)
call SIS_call_tracer_column_fns(dt_slow, G, IG, CS%tracer_flow_CSp, IST%mH_ice, mi_old)
call disable_SIS_averaging(CS%diag)
call accumulate_bottom_input(IST, OSS, FIA, IOF, dt_slow, G, IG, CS%sum_output_CSp)
if (CS%column_check) &
call write_ice_statistics(IST, CS%Time, CS%n_calls, G, IG, CS%sum_output_CSp, &
message=" Post_thermo A", check_column=.true.)
call adjust_ice_categories(IST%mH_ice, IST%mH_snow, IST%mH_pond, IST%part_size, &
IST%TrReg, G, IG, CS%ice_transport_CSp) !Niki: add ridging?
if (CS%column_check) &
call write_ice_statistics(IST, CS%Time, CS%n_calls, G, IG, CS%sum_output_CSp, &
message=" Post_thermo B ", check_column=.true.)
end subroutine slow_thermodynamics
subroutine SIS2_thermodynamics(IST, dt_slow, CS, OSS, FIA, IOF, G, IG)
type(ice_state_type), intent(inout) :: IST
real, intent(in) :: dt_slow ! The thermodynamic step, in s.
type(slow_thermo_CS), pointer :: CS
type(ocean_sfc_state_type), intent(in) :: OSS
type(fast_ice_avg_type), intent(inout) :: FIA
type(ice_ocean_flux_type), intent(inout) :: IOF
type(SIS_hor_grid_type), intent(inout) :: G
type(ice_grid_type), intent(inout) :: IG
! This subroutine does the thermodynamic calculations in the same order as SIS1,
! but with a greater emphasis on enthalpy as the dominant state variable.
real, dimension(G%isc:G%iec,G%jsc:G%jec) :: salt_change, h2o_change, bsnk, tmp2d
real, dimension(SZI_(G),SZJ_(G),1:IG%CatIce) :: snow_to_ice
real, dimension(G%isc:G%iec,G%jsc:G%jec) :: Obs_sst, Obs_h_ice ! for qflux calculation
real, dimension(G%isc:G%iec,G%jsc:G%jec,2) :: Obs_cn_ice ! partition 2 = ice concentration
real, dimension(G%isc:G%iec,G%jsc:G%jec) :: icec, icec_obs
real, dimension(SZI_(G),SZJ_(G)) :: &
qflx_lim_ice, qflx_res_ice, &
cool_nudge, & ! A heat flux out of the sea ice that
! acts to create sea-ice, in W m-2.
net_melt ! The net mass flux from the ice and snow into the
! ocean due to melting and freezing integrated
! across all categories, in kg m-2 s-1.
real, dimension(SZI_(G),SZJ_(G),1:IG%CatIce) :: heat_in, enth_prev, enth
real, dimension(SZI_(G),SZJ_(G)) :: heat_in_col, enth_prev_col, enth_col, enth_mass_in_col
real, dimension(IG%NkIce) :: S_col ! The salinity of a column of ice, in g/kg.
real, dimension(IG%NkIce+1) :: Salin ! The conserved bulk salinity of each
! layer in g/kg, with the salinity of
! newly formed ice in layer NkIce+1.
real, dimension(0:IG%NkIce) :: m_lay ! The masses of a column of ice and snow, in kg m-2.
real, dimension(0:IG%NkIce) :: Tcol0 ! The temperature of a column of ice and snow, in degC.
real, dimension(0:IG%NkIce) :: S_col0 ! The salinity of a column of ice and snow, in g/kg.
real, dimension(0:IG%NkIce) :: Tfr_col0 ! The freezing temperature of a column of ice and snow, in degC.
real, dimension(0:IG%NkIce+1) :: &
enthalpy ! The initial enthalpy of a column of ice and snow
! and the surface ocean, in enth_units (often J/kg).
real, dimension(IG%CatIce) :: frazil_cat ! The frazil heating applied to each thickness
! category, averaged over the area of that category in J m-2.
real :: enthalpy_ocean ! The enthalpy of the ocean surface waters, in Enth_units.
real :: heat_fill_val ! An enthalpy to use for massless categories, in enth_units.
real :: T_freeze_surf ! The freezing temperature at the surface salinity of
! the ocean, in deg C.
real :: I_part ! The inverse of a part_size, nondim.
logical :: spec_thermo_sal ! If true, use the specified salinities of the
! various sub-layers of the ice for all thermodynamic
! calculations; otherwise use the prognostic
! salinity fields for these calculations.
type(EOS_type), pointer :: EOS
real :: Cp_water
real :: drho_dT(1), drho_dS(1), pres_0(1)
real :: rho_ice ! The nominal density of sea ice in kg m-3.
real :: Idt_slow ! The inverse of the thermodynamic step, in s-1.
real :: yr_dtslow ! The ratio of 1 year to the thermodyamic time step, used
! to change the units of several diagnostics to rate yr-1
real :: heat_to_ocn, h2o_ice_to_ocn, h2o_ocn_to_ice, evap_from_ocn, sn2ic, bablt
real :: salt_to_ice ! The flux of salt from the ocean to the ice, in kg m-2 s-1.
! This may be of either sign; in some places it is an
! average over the whole cell, while in others just a partition.
real :: mtot_ice ! The total mass of ice ans snow in a cell, in kg m-2.
real :: e2m_tot ! The total enthalpy requred to melt all ice and snow, in J m-2.
real :: enth_evap, enth_ice_to_ocn, enth_ocn_to_ice, enth_snowfall
real :: tot_heat, heating, tot_frazil, heat_mass_in, heat_input
real :: mass_in, mass_here, mass_prev, mass_imb
real :: enth_units, I_enth_units ! The units of enthaply and their inverse.
real :: frac_keep, frac_melt ! The fraction of ice and snow to keep or remove, nd.
real :: ice_melt_lay ! The amount of excess ice removed from each layer in kg/m2.
real :: snow_melt ! The amount of excess snow that is melted, in kg/m2.
real :: enth_freeze ! The freezing point enthalpy of a layer, in enth_units.
real :: enth_to_melt ! The enthalpy addition required to melt the excess ice
! and snow in enth_unit kg/m2.
real :: I_Nk ! The inverse of the number of layers in the ice, nondim.
real :: kg_H_Nk ! The conversion factor from units of H to kg/m2 over Nk.
real :: part_sum ! A running sum of partition sizes.
real :: d_enth ! The change in enthalpy between categories.
real :: fill_frac ! The fraction of the difference between the thicknesses
! in thin categories that will be removed within a single
! timestep with filling_frazil.
integer :: i, j, k, l, m, n, isc, iec, jsc, jec, ncat, NkIce
integer :: k_merge
real :: LatHtFus ! The latent heat of fusion of ice in J/kg.
real :: LatHtVap ! The latent heat of vaporization of water at 0C in J/kg.
real, parameter :: T_0degC = 273.15 ! 0 degrees C in Kelvin
real :: tot_heat_in, enth_here, enth_imb, norm_enth_imb, emic2, tot_heat_in2, enth_imb2
isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; ncat = IG%CatIce
NkIce = IG%NkIce ; I_Nk = 1.0 / NkIce ; kg_H_Nk = IG%H_to_kg_m2 * I_Nk
Idt_slow = 0.0 ; if (dt_slow > 0.0) Idt_slow = 1.0/dt_slow
call get_SIS2_thermo_coefs(IST%ITV, ice_salinity=S_col, enthalpy_units=enth_units, &
rho_ice=rho_ice, specified_thermo_salinity=spec_thermo_sal, &
Latent_fusion=LatHtFus, Latent_vapor=LatHtVap)
S_col0(0) = 0.0 ; do m=1,NkIce ; S_col0(m) = S_col(m) ; enddo
call calculate_T_Freeze(S_col0(0:NkIce), Tfr_col0(0:NkIce), IST%ITV)
I_enth_units = 1.0 / enth_units
heat_fill_val = Enth_from_TS(0.0, 0.0, IST%ITV)
if (.not.spec_thermo_sal) call SIS_error(FATAL, "SIS2_thermodynamics is not "//&
"prepared for SPECIFIED_THERMO_SALINITY to be false.")
call mpp_clock_begin(iceClock6)
! Add any heat fluxes to restore the sea-ice properties toward a prescribed
! state, potentially including freshwater fluxes to avoid driving oceanic
! convection.
if (CS%nudge_sea_ice) then
if (.not.allocated(IOF%melt_nudge)) allocate(IOF%melt_nudge(isc:iec,jsc:jec))
cool_nudge(:,:) = 0.0 ; IOF%melt_nudge(:,:) = 0.0
icec(:,:) = 0.0
call data_override('ICE','icec',icec_obs,CS%Time)
do k=1,ncat ; do j=jsc,jec ; do i=isc,iec
icec(i,j) = icec(i,j) + IST%part_size(i,j,k)
enddo ; enddo ; enddo
pres_0(:) = 0.0
call get_SIS2_thermo_coefs(IST%ITV, Cp_SeaWater=Cp_water, EOS=EOS)
do j=jsc,jec ; do i=isc,iec
if (icec(i,j) < icec_obs(i,j) - CS%nudge_conc_tol) then
cool_nudge(i,j) = CS%nudge_sea_ice_rate * &
((icec_obs(i,j)-CS%nudge_conc_tol) - icec(i,j))**2.0 ! W/m2
if (CS%nudge_stab_fac /= 0.0) then
if (OSS%t_ocn(i,j) > T_Freeze(OSS%s_surf(i,j), IST%ITV) ) then
call calculate_density_derivs(OSS%t_ocn(i:i,j),OSS%s_surf(i:i,j),pres_0,&
drho_dT,drho_dS,1,1,EOS)
IOF%melt_nudge(i,j) = CS%nudge_stab_fac * (-cool_nudge(i,j)*drho_dT(1)) / &
((Cp_Water*drho_dS(1)) * max(OSS%s_surf(i,j), 1.0) )
endif
endif
elseif (icec(i,j) > icec_obs(i,j) + CS%nudge_conc_tol) then
! Heat the ice but do not apply a fresh water flux.
cool_nudge(i,j) = -CS%nudge_sea_ice_rate * &
(icec(i,j) - (icec_obs(i,j)+CS%nudge_conc_tol))**2.0 ! W/m2
endif
if (cool_nudge(i,J) > 0.0) then
FIA%frazil_left(i,j) = FIA%frazil_left(i,j) + cool_nudge(i,j)*dt_slow
elseif (cool_nudge(i,J) < 0.0) then
FIA%bheat(i,j) = FIA%bheat(i,j) - cool_nudge(i,j)
endif
enddo ; enddo
endif
if (CS%do_ice_restore) then
! get observed ice thickness for ice restoring, if calculating qflux
call get_sea_surface(CS%Time, Obs_SST, Obs_cn_ice, Obs_h_ice)
call get_SIS2_thermo_coefs(IST%ITV, Latent_fusion=LatHtFus)
qflx_res_ice(:,:) = 0.0
do j=jsc,jec ; do i=isc,iec
e2m_tot = 0.0
! Calculate the average ice and snow enthalpy relative to freezing per unit area.
do k=1,ncat
if (IST%part_size(i,j,k)*IST%mH_ice(i,j,k)>0.0) then
e2m_tot = (IST%part_size(i,j,k)*IST%mH_snow(i,j,k)) * IG%H_to_kg_m2 * &
((enthalpy_liquid_freeze(0.0, IST%ITV) - &
IST%enth_snow(i,j,k,1)) * I_enth_units)
if (spec_thermo_sal) then ; do m=1,NkIce
e2m_tot = e2m_tot + (IST%part_size(i,j,k)*IST%mH_ice(i,j,k)*kg_H_Nk) * &
((enthalpy_liquid_freeze(S_col(m), IST%ITV) - &
IST%enth_ice(i,j,k,m)) * I_enth_units)
enddo ; else ; do m=1,NkIce
e2m_tot = e2m_tot + (IST%part_size(i,j,k)*IST%mH_ice(i,j,k)*kg_H_Nk) * &
((enthalpy_liquid_freeze(IST%sal_ice(i,j,k,m), IST%ITV) - &
IST%enth_ice(i,j,k,m)) * I_enth_units)
enddo ; endif
endif
enddo
qflx_res_ice(i,j) = -(LatHtFus*Rho_ice*Obs_h_ice(i,j)*Obs_cn_ice(i,j,2)-e2m_tot) / &
(86400.0*CS%ice_restore_timescale)
if (qflx_res_ice(i,j) < 0.0) then
FIA%frazil_left(i,j) = FIA%frazil_left(i,j) - qflx_res_ice(i,j)*dt_slow
elseif (qflx_res_ice(i,j) > 0.0) then
FIA%bheat(i,j) = FIA%bheat(i,j) + qflx_res_ice(i,j)
endif
enddo ; enddo
endif
call mpp_clock_end(iceClock6)
if (CS%column_check) then
enth_prev(:,:,:) = 0.0 ; heat_in(:,:,:) = 0.0
enth_prev_col(:,:) = 0.0 ; heat_in_col(:,:) = 0.0 ; enth_mass_in_col(:,:) = 0.0
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,ncat,IST,enth_prev,G,I_Nk, &
!$OMP heat_in_col,dt_slow,enth_prev_col,NkIce,FIA)
do j=jsc,jec
do k=1,ncat ; do i=isc,iec ; if (IST%mH_ice(i,j,k) > 0.0) then
enth_prev(i,j,k) = IST%mH_snow(i,j,k) * IST%enth_snow(i,j,k,1)
do m=1,NkIce
enth_prev(i,j,k) = enth_prev(i,j,k) + (IST%mH_ice(i,j,k)*I_Nk) * IST%enth_ice(i,j,k,m)
enddo
endif ; enddo ; enddo
do i=isc,iec
heat_in_col(i,j) = heat_in_col(i,j) - FIA%frazil_left(i,j)
heat_in_col(i,j) = heat_in_col(i,j) - IST%part_size(i,j,0) * dt_slow*FIA%flux_t_top(i,j,0)
enddo
do k=1,ncat ; do i=isc,iec
if (IST%part_size(i,j,k) > 0.0) then
heat_in_col(i,j) = heat_in_col(i,j) + IST%part_size(i,j,k) * &
(FIA%tmelt(i,j,k) + FIA%bmelt(i,j,k) - dt_slow*FIA%bheat(i,j))
endif
if (IST%part_size(i,j,k)*IST%mH_ice(i,j,k) > 0.0) then
enth_prev_col(i,j) = enth_prev_col(i,j) + &
(IST%mH_snow(i,j,k)*IST%part_size(i,j,k)) * IST%enth_snow(i,j,k,1)
do m=1,NkIce
enth_prev_col(i,j) = enth_prev_col(i,j) + &
(IST%mH_ice(i,j,k)*IST%part_size(i,j,k)*I_Nk) * IST%enth_ice(i,j,k,m)
enddo
endif
enddo ; enddo
enddo
endif
call mpp_clock_begin(iceClock5)
snow_to_ice(:,:,:) = 0.0 ; net_melt(:,:) = 0.0
bsnk(:,:) = 0.0
salt_change(:,:) = 0.0
h2o_change(:,:) = 0.0
!$OMP parallel default(none) shared(isc,iec,jsc,jec,ncat,G,IST,salt_change,kg_H_Nk, &
!$OMP h2o_change,NkIce,IG,CS,IOF,FIA)
if (CS%ice_rel_salin <= 0.0) then
!$OMP do
do j=jsc,jec ; do m=1,NkIce ; do k=1,ncat ; do i=isc,iec
salt_change(i,j) = salt_change(i,j) - &
(IST%sal_ice(i,j,k,m)*(IST%mH_ice(i,j,k)*kg_H_Nk)) * IST%part_size(i,j,k)
enddo ; enddo ; enddo ; enddo
endif
!$OMP do
do j=jsc,jec ; do k=1,ncat ; do i=isc,iec
h2o_change(i,j) = h2o_change(i,j) - IST%part_size(i,j,k) * &
IG%H_to_kg_m2*(IST%mH_snow(i,j,k) + IST%mH_ice(i,j,k))
enddo ; enddo ; enddo
! Start accumulating the fluxes at the ocean's surface.
!$OMP do
do j=jsc,jec ; do i=isc,iec
IOF%flux_t_ocn_top(i,j) = IST%part_size(i,j,0) * FIA%flux_t_top(i,j,0)
IOF%flux_q_ocn_top(i,j) = IST%part_size(i,j,0) * FIA%flux_q_top(i,j,0)
IOF%flux_lw_ocn_top(i,j) = IST%part_size(i,j,0) * FIA%flux_lw_top(i,j,0)
IOF%flux_lh_ocn_top(i,j) = IST%part_size(i,j,0) * FIA%flux_lh_top(i,j,0)
IOF%flux_sw_vis_dir_ocn(i,j) = IST%part_size(i,j,0) * FIA%flux_sw_vis_dir_top(i,j,0)
IOF%flux_sw_vis_dif_ocn(i,j) = IST%part_size(i,j,0) * FIA%flux_sw_vis_dif_top(i,j,0)
IOF%flux_sw_nir_dir_ocn(i,j) = IST%part_size(i,j,0) * FIA%flux_sw_nir_dir_top(i,j,0)
IOF%flux_sw_nir_dif_ocn(i,j) = IST%part_size(i,j,0) * FIA%flux_sw_nir_dif_top(i,j,0)
IOF%lprec_ocn_top(i,j) = IST%part_size(i,j,0) * FIA%lprec_top(i,j,0)
IOF%fprec_ocn_top(i,j) = IST%part_size(i,j,0) * FIA%fprec_top(i,j,0)
enddo ; enddo
! mw/new precip will eventually be intercepted by pond eliminating need for next 3 lines
!$OMP do
do j=jsc,jec ; do k=1,ncat ; do i=isc,iec
IOF%lprec_ocn_top(i,j) = IOF%lprec_ocn_top(i,j) + &
IST%part_size(i,j,k) * FIA%lprec_top(i,j,k)
enddo ; enddo ; enddo
!$OMP end parallel
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,ncat,G,IST,S_col0,NkIce,S_col, &
!$OMP dt_slow,snow_to_ice,heat_in,I_NK,enth_units, &
!$OMP enth_prev,enth_mass_in_col,Idt_slow,bsnk, &
!$OMP salt_change,net_melt,kg_H_nk,LatHtFus,LatHtVap,&
!$OMP IG,CS,OSS,FIA,IOF) &
!$OMP private(mass_prev,enthalpy,enthalpy_ocean,Salin, &
!$OMP heat_to_ocn,h2o_ice_to_ocn,h2o_ocn_to_ice, &
!$OMP evap_from_ocn,salt_to_ice,bablt,enth_evap, &
!$OMP enth_ice_to_ocn,enth_ocn_to_ice,heat_input, &
!$OMP heat_mass_in,mass_in,mass_here,enth_here, &
!$OMP tot_heat_in,enth_imb,mass_imb,norm_enth_imb, &
!$OMP m_lay, mtot_ice, &
!$OMP T_Freeze_surf,I_part,sn2ic,enth_snowfall)
do j=jsc,jec ; do k=1,ncat ; do i=isc,iec
if (G%mask2dT(i,j) > 0 .and. IST%part_size(i,j,k) > 0) then
! reshape the ice based on fluxes
if (CS%column_check) then
mass_prev = IST%mH_snow(i,j,k)
mass_prev = mass_prev + IST%mH_ice(i,j,k)
endif
! evap_from_ocn = 0.0; h2o_ice_to_ocn = 0.0; heat_to_ocn = 0.0
if (IST%mH_snow(i,j,k) == 0.0) IST%enth_snow(i,j,k,1) = &
enth_from_TS(Temp_from_En_S(IST%enth_ice(i,j,k,1), S_col0(1), IST%ITV), &
0.0, IST%ITV)
enthalpy(0) = IST%enth_snow(i,j,k,1)
do m=1,NkIce ; enthalpy(m) = IST%enth_ice(i,j,k,m) ; enddo
enthalpy_ocean = enthalpy_liquid(OSS%t_ocn(i,j), OSS%s_surf(i,j), IST%ITV)
enthalpy(NkIce+1) = enthalpy_ocean
if (CS%ice_rel_salin > 0.0) then
do m=1,NkIce ; Salin(m) = IST%sal_ice(i,j,k,m) ; enddo
salin(NkIce+1) = CS%ice_rel_salin * OSS%s_surf(i,j)
else
do m=1,NkIce+1 ; Salin(m) = CS%ice_bulk_salin ; enddo
endif
m_lay(0) = IST%mH_snow(i,j,k) * IG%H_to_kg_m2
do m=1,NkIce ; m_lay(m) = IST%mH_ice(i,j,k) * kg_H_Nk ; enddo
! mw/new - melt pond size is now adjusted here (rain ignored in resize, for now)
call ice_resize_SIS2(1-IST%part_size(i,j,0), IST%mH_pond(i,j,k), m_lay, &
enthalpy, S_col, Salin, FIA%fprec_top(i,j,k)*dt_slow, &
FIA%lprec_top(i,j,k)*dt_slow, FIA%flux_q_top(i,j,k)*dt_slow, &
FIA%tmelt(i,j,k), FIA%bmelt(i,j,k), NkIce, &
heat_to_ocn, h2o_ice_to_ocn, h2o_ocn_to_ice, evap_from_ocn, &
snow_to_ice(i,j,k), salt_to_ice, IST%ITV, CS%ice_thm_CSp, bablt, &
enth_evap, enth_ice_to_ocn, enth_ocn_to_ice)
IST%mH_snow(i,j,k) = m_lay(0) * IG%kg_m2_to_H
call rebalance_ice_layers(m_lay, mtot_ice, Enthalpy, Salin, NkIce)
IST%mH_ice(i,j,k) = mtot_ice * IG%kg_m2_to_H
if (IST%mH_ice(i,j,k) == 0.0) then
do m=1,NkIce ; IST%enth_ice(i,j,k,m) = enthalpy_liquid_freeze(S_col(m), IST%ITV) ; enddo
else
do m=1,NkIce ; IST%enth_ice(i,j,k,m) = enthalpy(m) ; enddo
endif
if (CS%ice_rel_salin > 0.0) then
if (IST%mH_ice(i,j,k) == 0.0) then
do m=1,NkIce ; IST%sal_ice(i,j,k,m) = 0.0 ; enddo
else
do m=1,NkIce ; IST%sal_ice(i,j,k,m) = Salin(m) ; enddo
endif
salt_change(i,j) = salt_change(i,j) + IST%part_size(i,j,k) * salt_to_ice
endif
! The snow enthalpy should not have changed. This should do nothing.
! IST%enth_snow(i,j,k,1) = Enthalpy(0)
enth_snowfall = ((dt_slow*FIA%fprec_top(i,j,k)) * enthalpy(0))
IOF%Enth_Mass_in_atm(i,j) = IOF%Enth_Mass_in_atm(i,j) + &
IST%part_size(i,j,k) * enth_snowfall
! IOF%Enth_Mass_in_ocn(i,j) = IOF%Enth_Mass_in_ocn(i,j) + &
! IST%part_size(i,j,k) * enth_ocn_to_ice
IOF%Enth_Mass_in_ocn(i,j) = IOF%Enth_Mass_in_ocn(i,j) + &
IST%part_size(i,j,k) * (h2o_ocn_to_ice * enthalpy_ocean)
IOF%Enth_Mass_out_ocn(i,j) = IOF%Enth_Mass_out_ocn(i,j) - &
IST%part_size(i,j,k) * enth_ice_to_ocn
IOF%Enth_Mass_out_atm(i,j) = IOF%Enth_Mass_out_atm(i,j) - &
IST%part_size(i,j,k) * enth_evap
if (CS%column_check) then
heat_in(i,j,k) = heat_in(i,j,k) + FIA%tmelt(i,j,k) + FIA%bmelt(i,j,k) - &
(heat_to_ocn - (LatHtVap+LatHtFus)*evap_from_ocn)
heat_input = FIA%tmelt(i,j,k) + FIA%bmelt(i,j,k) - (heat_to_ocn - (LatHtVap+LatHtFus)*evap_from_ocn)
heat_mass_in = enth_snowfall + enth_ocn_to_ice - enth_ice_to_ocn - enth_evap
mass_in = dt_slow*FIA%fprec_top(i,j,k) & ! +FIA%lprec_top(i,j,k) <- eventually
+ h2o_ocn_to_ice - h2o_ice_to_ocn &
- (dt_slow*FIA%flux_q_top(i,j,k)-evap_from_ocn)
mass_here = IST%mH_snow(i,j,k) + IST%mH_pond(i,j,k) + IST%mH_ice(i,j,k)
enth_here = IST%mH_snow(i,j,k) * IST%enth_snow(i,j,k,1)
do m=1,NkIce
enth_here = enth_here + (IST%mH_ice(i,j,k)*I_Nk) * IST%enth_ice(i,j,k,m)
enddo
tot_heat_in = IG%kg_m2_to_H*(enth_units*heat_input + heat_mass_in)
mass_in = mass_in*IG%kg_m2_to_H
enth_imb = enth_here - (enth_prev(i,j,k) + tot_heat_in)
mass_imb = mass_here - (mass_prev + mass_in)
if (abs(enth_imb) > CS%imb_tol * &
(abs(enth_here) + abs(enth_prev(i,j,k)) + abs(tot_heat_in)) ) then
norm_enth_imb = enth_imb / (abs(enth_here) + abs(enth_prev(i,j,k)) + abs(tot_heat_in))
enth_imb = enth_here - (enth_prev(i,j,k) + tot_heat_in)
endif
enth_mass_in_col(i,j) = enth_mass_in_col(i,j) + IST%part_size(i,j,k) * &
(enth_snowfall + enth_ocn_to_ice - enth_ice_to_ocn - enth_evap)
endif
IOF%flux_q_ocn_top(i,j) = IOF%flux_q_ocn_top(i,j) + IST%part_size(i,j,k) * &
(evap_from_ocn*Idt_slow) ! no ice, evaporation left
IOF%flux_lh_ocn_top(i,j) = IOF%flux_lh_ocn_top(i,j) + IST%part_size(i,j,k) * &
((LatHtVap*evap_from_ocn)*Idt_slow)
IOF%flux_t_ocn_top(i,j) = IOF%flux_t_ocn_top(i,j) + IST%part_size(i,j,k) * &
(FIA%bheat(i,j) - (heat_to_ocn - LatHtFus*evap_from_ocn)*Idt_slow)
IOF%flux_sw_vis_dif_ocn(i,j) = IOF%flux_sw_vis_dif_ocn(i,j) + IST%part_size(i,j,k) * &
(((FIA%flux_sw_vis_dir_top(i,j,k) + FIA%flux_sw_vis_dif_top(i,j,k)) + &
(FIA%flux_sw_nir_dir_top(i,j,k) + FIA%flux_sw_nir_dif_top(i,j,k))) * &
FIA%sw_abs_ocn(i,j,k))
net_melt(i,j) = net_melt(i,j) + IST%part_size(i,j,k) * &
((h2o_ice_to_ocn-h2o_ocn_to_ice)*Idt_slow)
bsnk(i,j) = bsnk(i,j) - IST%part_size(i,j,k)*bablt ! bot. melt. ablation
endif ! Applying surface fluxes to each category.
enddo ; enddo ; enddo
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,ncat,G,IG,IST,S_col0,NkIce,S_col, &
!$OMP dt_slow,snow_to_ice,heat_in,I_NK,enth_units, &
!$OMP enth_prev,enth_mass_in_col,Idt_slow,bsnk, &
!$OMP salt_change,net_melt,kg_h_Nk,LatHtFus,FIA,CS,OSS,IOF) &
!$OMP private(mass_prev,enthalpy,enthalpy_ocean,Salin, &
!$OMP heat_to_ocn,h2o_ice_to_ocn,h2o_ocn_to_ice, &
!$OMP evap_from_ocn,salt_to_ice,bablt,enth_evap, &
!$OMP enth_ice_to_ocn,enth_ocn_to_ice,heat_input, &
!$OMP heat_mass_in,mass_in,mass_here,enth_here, &
!$OMP tot_heat_in,enth_imb,mass_imb,norm_enth_imb, &
!$OMP m_lay,mtot_ice,frazil_cat,k_merge,part_sum, &
!$OMP fill_frac,d_enth, &
!$OMP T_Freeze_surf,I_part,sn2ic,enth_snowfall)
do j=jsc,jec ; do i=isc,iec ; if (FIA%frazil_left(i,j)>0.0) then
frazil_cat(1:ncat) = 0.0
k_merge = 1 ! Find the category that will be combined with the ice free category.
if (.not.CS%filling_frazil) then
do k=1,ncat ; if (IST%part_size(i,j,0)+IST%part_size(i,j,k)>0.01) then
! absorb frazil in thinest ice partition available (was ...>0.0)
! raised above threshold from 0 to 0.01 to avert ocean-ice model blow-ups
k_merge = k ; exit
endif ; enddo
endif
! if (IST%part_size(i,j,0) > 0.0) then
k = k_merge
T_Freeze_surf = T_Freeze(OSS%s_surf(i,j),IST%ITV)
! Combine the ice-free part size with one of the categories.
I_part = 1.0 / (IST%part_size(i,j,k) + IST%part_size(i,j,0))
IST%mH_snow(i,j,k) = (IST%mH_snow(i,j,k) * IST%part_size(i,j,k)) * I_part
IST%mH_ice(i,j,k) = (IST%mH_ice(i,j,k) * IST%part_size(i,j,k)) * I_part
IST%t_surf(i,j,k) = (IST%t_surf(i,j,k) * IST%part_size(i,j,k) + &
(T_0degC + T_Freeze_surf)*IST%part_size(i,j,0)) * I_part
IST%part_size(i,j,k) = IST%part_size(i,j,k) + IST%part_size(i,j,0)
IST%part_size(i,j,0) = 0.0
! endif
if (CS%filling_frazil) then
if (CS%fraz_fill_time < 0.0) then
frazil_cat(k) = FIA%frazil_left(i,J)
FIA%frazil_left(i,j) = 0.0
else
part_sum = 0.0
fill_frac = 1.0 ; if (CS%fraz_fill_time > 0.0) &
fill_frac = dt_slow / (dt_slow + CS%fraz_fill_time)
do k=1,ncat-1
part_sum = part_sum + IST%part_size(i,j,k)
d_enth = fill_frac * max(0.0, LatHtFus * IG%H_to_kg_m2 * &
(IG%mH_cat_bound(k+1) - IST%mH_ice(i,j,k)))
if (d_enth*part_sum > FIA%frazil_left(i,j)) then
frazil_cat(k) = FIA%frazil_left(i,j) / part_sum
FIA%frazil_left(i,j) = 0.0
exit
else
frazil_cat(k) = d_enth
FIA%frazil_left(i,j) = FIA%frazil_left(i,j) - frazil_cat(k)*part_sum
endif
enddo
if (FIA%frazil_left(i,j) > 0.0) &
frazil_cat(ncat) = FIA%frazil_left(i,j)
! Note that at this point we should have that part_sum = 1.0.
do k=ncat-1,1 ; frazil_cat(k) = frazil_cat(k) + frazil_cat(k+1) ; enddo
endif
else
! Set the frazil that is absorbed in this category and remove it from
! the overall frazil energy.
I_part = 1.0 / (IST%part_size(i,j,k))
frazil_cat(k_merge) = FIA%frazil_left(i,j) * I_part
FIA%frazil_left(i,j) = 0.0
endif
do k=1,ncat ; if (frazil_cat(k) > 0.0) then
if (CS%column_check) then
enth_prev(i,j,k) = 0.0 ; heat_in(i,j,k) = 0.0
if (IST%mH_ice(i,j,k) > 0.0) then
enth_prev(i,j,k) = IST%mH_snow(i,j,k) * IST%enth_snow(i,j,k,1)
do m=1,NkIce
enth_prev(i,j,k) = enth_prev(i,j,k) + (IST%mH_ice(i,j,k)*I_Nk) * IST%enth_ice(i,j,k,m)
enddo
enth_prev(i,j,k) = enth_prev(i,j,k) * IST%part_size(i,j,k)
endif
endif
! Set up the columns of enthalpy, salinity, and mass.
enthalpy(0) = IST%enth_snow(i,j,k,1)
do m=1,NkIce ; enthalpy(m) = IST%enth_ice(i,j,k,m) ; enddo
enthalpy_ocean = enthalpy_liquid(OSS%t_ocn(i,j), OSS%s_surf(i,j), IST%ITV)
enthalpy(NkIce+1) = enthalpy_ocean
if (CS%ice_rel_salin > 0.0) then
do m=1,NkIce ; Salin(m) = IST%sal_ice(i,j,k,m) ; enddo
salin(NkIce+1) = CS%ice_rel_salin * OSS%s_surf(i,j)
else
do m=1,NkIce+1 ; Salin(m) = CS%ice_bulk_salin ; enddo
endif
m_lay(0) = IST%mH_snow(i,j,k) * IG%H_to_kg_m2
do m=1,NkIce ; m_lay(m) = IST%mH_ice(i,j,k) * kg_H_Nk ; enddo
call add_frazil_SIS2(m_lay, enthalpy, S_col, Salin, frazil_cat(k), &
T_Freeze_surf, NkIce, h2o_ocn_to_ice, salt_to_ice, IST%ITV, &
CS%ice_thm_CSp, Enthalpy_freeze=enth_ocn_to_ice)
call rebalance_ice_layers(m_lay, mtot_ice, Enthalpy, Salin, NkIce)
! Unpack the columns of mass, enthalpy and salinity.
IST%mH_snow(i,j,k) = m_lay(0) * IG%kg_m2_to_H
IST%mH_ice(i,j,k) = mtot_ice * IG%kg_m2_to_H
if (IST%mH_ice(i,j,k) == 0.0) then
do m=1,NkIce ; IST%enth_ice(i,j,k,m) = enthalpy_liquid_freeze(S_col(m), IST%ITV) ; enddo
else
do m=1,NkIce ; IST%enth_ice(i,j,k,m) = enthalpy(m) ; enddo
endif
if (CS%ice_rel_salin > 0.0) then
if (IST%mH_ice(i,j,k) == 0.0) then
do m=1,NkIce ; IST%sal_ice(i,j,k,m) = 0.0 ; enddo
else
do m=1,NkIce ; IST%sal_ice(i,j,k,m) = Salin(m) ; enddo
endif
salt_change(i,j) = salt_change(i,j) + IST%part_size(i,j,k) * salt_to_ice
endif
! IOF%Enth_Mass_in_ocn(i,j) = IOF%Enth_Mass_in_ocn(i,j) + &
! IST%part_size(i,j,k) * enth_ocn_to_ice
IOF%Enth_Mass_in_ocn(i,j) = IOF%Enth_Mass_in_ocn(i,j) + &
IST%part_size(i,j,k) * (h2o_ocn_to_ice * enthalpy_ocean)
net_melt(i,j) = net_melt(i,j) - &
(h2o_ocn_to_ice * IST%part_size(i,j,k)) * Idt_slow
if (CS%column_check) then
heat_in(i,j,k) = heat_in(i,j,k) - frazil_cat(k)
enth_here = IST%mH_snow(i,j,k) * IST%enth_snow(i,j,k,1)
do m=1,NkIce
enth_here = enth_here + (IST%mH_ice(i,j,k)*I_Nk) * IST%enth_ice(i,j,k,m)
enddo
enth_here = enth_here * IST%part_size(i,j,k)
tot_heat_in = (enth_units * heat_in(i,j,k) + enth_ocn_to_ice) * IST%part_size(i,j,k)
enth_imb = enth_here - (enth_prev(i,j,k) + tot_heat_in)
if (abs(enth_imb) > CS%imb_tol * &
(abs(enth_here) + abs(enth_prev(i,j,k)) + abs(tot_heat_in)) ) then
enth_imb = enth_here - (enth_prev(i,j,k) + tot_heat_in)
endif
enth_mass_in_col(i,j) = enth_mass_in_col(i,j) + IST%part_size(i,j,k) * enth_ocn_to_ice
endif
endif ; enddo ! frazil_cat>0, k-loop
endif ; enddo ; enddo ! frazil>0, i-, and j-loops