-
Notifications
You must be signed in to change notification settings - Fork 45
/
Copy pathice_model.F90
2835 lines (2407 loc) · 142 KB
/
ice_model.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
!> This is the central module for the SIS2 sea ice model.
module ice_model_mod
! This file is a part of SIS2. See LICENSE.md for the license.
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
! 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 manages fluxes between sub-modules, many diagnostics, and the !
! overall time stepping of the sea ice. Sea ice dynamics are handled in !
! ice_dyn_bgrid.F90 or ice_dyn_cgrid.F90, while the transport of mass, heat, !
! and tracers occurs in ice_transport.F90. Sea ice thermodynamics is treated !
! in ice_thm.F90 and other modules that are subsequently called from there. !
! The Lagrangian icebergs code of Adcroft and Martin is called from SIS. !
! The original SIS was developed by Mike Winton (Michael.Winton@noaa.gov). !
! SIS2 has been developed by Robert Hallberg and Mike Winton, with !
! contributions from many people at NOAA/GFDL, including Alistair Adcroft and !
! Niki Zadeh. !
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
use MOM_cpu_clock, only : CLOCK_COMPONENT, CLOCK_SUBCOMPONENT
use MOM_domains, only : MOM_domain_type
use MOM_domains, only : pass_var, pass_vector, AGRID, BGRID_NE, CGRID_NE
use MOM_domains, only : fill_symmetric_edges, MOM_domains_init, clone_MOM_domain
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, read_param, param_file_type
use MOM_file_parser, only : open_param_file, close_param_file
use MOM_hor_index, only : hor_index_type, hor_index_init
use MOM_io, only : file_exists
use MOM_obsolete_params, only : obsolete_logical, obsolete_real
use MOM_string_functions, only : uppercase, extract_real
use MOM_time_manager, only : time_type, time_type_to_real, real_to_time
use MOM_time_manager, only : operator(+), operator(-)
use MOM_time_manager, only : operator(>), operator(*), operator(/), operator(/=)
use MOM_unit_scaling, only : unit_scale_type, unit_scaling_init, unit_scaling_end
use astronomy_mod, only : astronomy_init, astronomy_end
use astronomy_mod, only : universal_time, orbital_time, diurnal_solar, daily_mean_solar
use ocean_albedo_mod, only : compute_ocean_albedo ! ice sets ocean surface
use ocean_rough_mod, only : compute_ocean_roughness ! properties over water
use ice_bergs, only : icebergs, icebergs_run, icebergs_init, icebergs_end
use ice_boundary_types, only : ocean_ice_boundary_type, atmos_ice_boundary_type, land_ice_boundary_type
use ice_boundary_types, only : ocn_ice_bnd_type_chksum, atm_ice_bnd_type_chksum
use ice_boundary_types, only : lnd_ice_bnd_type_chksum
use ice_grid, only : set_ice_grid, ice_grid_end, ice_grid_type
use ice_spec_mod, only : get_sea_surface
use ice_type_mod, only : ice_data_type, dealloc_ice_arrays
use ice_type_mod, only : ice_type_slow_reg_restarts, ice_type_fast_reg_restarts
use ice_type_mod, only : Ice_public_type_chksum, Ice_public_type_bounds_check
use ice_type_mod, only : ice_model_restart, ice_stock_pe, ice_data_type_chksum
use SIS_ctrl_types, only : SIS_slow_CS, SIS_fast_CS
use SIS_ctrl_types, only : ice_diagnostics_init, ice_diags_fast_init
use SIS_debugging, only : chksum, uvchksum, Bchksum, SIS_debugging_init
use SIS_diag_mediator, only : set_SIS_axes_info, SIS_diag_mediator_init, SIS_diag_mediator_end
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_dyn_trans, only : SIS_dynamics_trans, SIS_multi_dyn_trans, update_icebergs
use SIS_dyn_trans, only : slab_ice_dyn_trans
use SIS_dyn_trans, only : SIS_dyn_trans_register_restarts, SIS_dyn_trans_init, SIS_dyn_trans_end
use SIS_dyn_trans, only : SIS_dyn_trans_read_alt_restarts, stresses_to_stress_mag
use SIS_dyn_trans, only : SIS_dyn_trans_transport_CS, SIS_dyn_trans_sum_output_CS
use SIS_fast_thermo, only : accumulate_deposition_fluxes, convert_frost_to_snow
use SIS_fast_thermo, only : do_update_ice_model_fast, avg_top_quantities, total_top_quantities
use SIS_fast_thermo, only : redo_update_ice_model_fast, find_excess_fluxes
use SIS_fast_thermo, only : infill_array, SIS_fast_thermo_init, SIS_fast_thermo_end
use SIS_framework, only : set_domain, nullify_domain, broadcast_domain
use SIS_restart, only : restore_SIS_state, query_initialized=>query_inited, SIS_restart_init
use SIS_restart, only : determine_is_new_run, is_new_run
use SIS_framework, only : coupler_1d_bc_type, coupler_2d_bc_type, coupler_3d_bc_type
use SIS_framework, only : coupler_type_spawn, coupler_type_initialized
use SIS_framework, only : coupler_type_rescale_data, coupler_type_copy_data
use SIS_fixed_initialization, only : SIS_initialize_fixed
use SIS_get_input, only : Get_SIS_input, directories
use SIS_hor_grid, only : SIS_hor_grid_type, set_hor_grid, SIS_hor_grid_end, set_first_direction
use SIS_open_boundary, only : ice_OBC_type
use SIS_optics, only : ice_optics_SIS2, SIS_optics_init, SIS_optics_end, SIS_optics_CS
use SIS_optics, only : VIS_DIR, VIS_DIF, NIR_DIR, NIR_DIF
use SIS_slow_thermo, only : slow_thermodynamics, SIS_slow_thermo_init, SIS_slow_thermo_end
use SIS_slow_thermo, only : SIS_slow_thermo_set_ptrs
use SIS_state_initialization, only : read_archaic_thermo_restarts, initialize_ice_categories
use SIS_state_initialization, only : ice_state_mass_init, ice_state_thermo_init
use SIS_sum_output, only : SIS_sum_output_init, write_ice_statistics
use SIS_tracer_flow_control, only : SIS_call_tracer_register, SIS_tracer_flow_control_init
use SIS_tracer_flow_control, only : SIS_tracer_flow_control_end
use SIS_tracer_registry, only : register_SIS_tracer, register_SIS_tracer_pair
use SIS_transcribe_grid, only : copy_dyngrid_to_SIS_horgrid, copy_SIS_horgrid_to_dyngrid
use SIS_transport, only : adjust_ice_categories
use SIS_types, only : ice_ocean_flux_type, alloc_ice_ocean_flux, dealloc_ice_ocean_flux
use SIS_types, only : ocean_sfc_state_type, alloc_ocean_sfc_state, dealloc_ocean_sfc_state, OSS_chksum
use SIS_types, only : fast_ice_avg_type, alloc_fast_ice_avg, dealloc_fast_ice_avg
use SIS_types, only : total_sfc_flux_type, alloc_total_sfc_flux, dealloc_total_sfc_flux
use SIS_types, only : ice_rad_type, ice_rad_register_restarts, dealloc_ice_rad, alloc_ice_rad
use SIS_types, only : simple_OSS_type, alloc_simple_OSS, dealloc_simple_OSS
use SIS_types, only : ice_state_type, alloc_IST_arrays, dealloc_IST_arrays
use SIS_types, only : IST_chksum, IST_bounds_check, ice_state_register_restarts
use SIS_types, only : ice_state_read_alt_restarts, register_fast_to_slow_restarts
use SIS_types, only : rescale_ice_state_restart_fields
use SIS_types, only : copy_IST_to_IST, copy_FIA_to_FIA, copy_sOSS_to_sOSS
use SIS_types, only : copy_TSF_to_TSF, redistribute_TSF_to_TSF, TSF_chksum
use SIS_types, only : copy_Rad_to_Rad, redistribute_Rad_to_Rad
use SIS_types, only : redistribute_IST_to_IST, redistribute_FIA_to_FIA
use SIS_types, only : redistribute_sOSS_to_sOSS, FIA_chksum, IOF_chksum, translate_OSS_to_sOSS
use SIS_utils, only : post_avg, ice_grid_chksum
use SIS2_ice_thm, only : ice_temp_SIS2, SIS2_ice_thm_init, SIS2_ice_thm_end
use SIS2_ice_thm, only : ice_thermo_init, ice_thermo_end, T_freeze, ice_thermo_type
use specified_ice, only : specified_ice_dynamics, specified_ice_init, specified_ice_CS
use specified_ice, only : specified_ice_end, specified_ice_sum_output_CS
implicit none ; private
#include <SIS2_memory.h>
public :: ice_data_type, ocean_ice_boundary_type, atmos_ice_boundary_type, land_ice_boundary_type
public :: ice_model_init, share_ice_domains, ice_model_end, ice_stock_pe
public :: update_ice_model_fast
public :: ice_model_restart ! for intermediate restarts
public :: ocn_ice_bnd_type_chksum, atm_ice_bnd_type_chksum
public :: lnd_ice_bnd_type_chksum, ice_data_type_chksum
public :: update_ice_atm_deposition_flux
public :: unpack_ocean_ice_boundary, unpack_ocn_ice_bdry, exchange_slow_to_fast_ice, set_ice_surface_fields
public :: unpack_ocean_ice_boundary_calved_shelf_bergs
public :: ice_model_fast_cleanup, unpack_land_ice_boundary
public :: exchange_fast_to_slow_ice, update_ice_model_slow
public :: update_ice_slow_thermo, update_ice_dynamics_trans
!>@{ CPU time clock IDs
integer :: iceClock
integer :: ice_clock_slow, ice_clock_fast, ice_clock_exchange
!!@}
integer, parameter :: REDIST=2 !< Redistribute for exchange
integer, parameter :: DIRECT=3 !< Use direct exchange
contains
!-----------------------------------------------------------------------
!> Update the sea-ice state due to slow processes, including dynamics,
!! freezing and melting, precipitation, and transport.
subroutine update_ice_model_slow(Ice)
type(ice_data_type), intent(inout) :: Ice !< The publicly visible ice data type.
call update_ice_slow_thermo(Ice)
call update_ice_dynamics_trans(Ice)
end subroutine update_ice_model_slow
!-----------------------------------------------------------------------
!> Update the sea-ice state due to slow thermodynamic processes, including
!! freezing and melting, precipitation, and brine drainage, and possibly also
!! also the accumulated effects of faster thermodynamic processes
subroutine update_ice_slow_thermo(Ice)
type(ice_data_type), intent(inout) :: Ice !< The publicly visible ice data type.
! These pointers are used to simplify the code below.
type(ice_grid_type), pointer :: sIG => NULL()
type(SIS_hor_grid_type), pointer :: sG => NULL()
type(ice_state_type), pointer :: sIST => NULL()
type(fast_ice_avg_type), pointer :: FIA => NULL()
type(ice_rad_type), pointer :: Rad => NULL()
type(unit_scale_type), pointer :: US => NULL()
real :: dt_slow ! The time step over which to advance the model [T ~> s].
integer :: i, j, i2, j2, i_off, j_off
if (.not.associated(Ice%sCS)) call SIS_error(FATAL, &
"The pointer to Ice%sCS must be associated in update_ice_slow_thermo.")
sIST => Ice%sCS%IST ; sIG => Ice%sCS%IG ; sG => Ice%sCS%G ; FIA => Ice%sCS%FIA
Rad => Ice%sCS%Rad ; US => Ice%sCS%US
call cpu_clock_begin(iceClock) ; call cpu_clock_begin(ice_clock_slow)
! Advance the slow PE clock to give the end time of the slow timestep. There
! is a separate clock inside the fCS that is advanced elsewhere.
Ice%sCS%Time = Ice%sCS%Time + Ice%sCS%Time_step_slow
if (.not.associated(Ice%fCS)) then
Ice%Time = Ice%sCS%Time
endif
dt_slow = US%s_to_T*time_type_to_real(Ice%sCS%Time_step_slow)
if (Ice%sCS%debug) then
call Ice_public_type_chksum("Start update_ice_slow_thermo", Ice, check_slow=.true.)
call FIA_chksum("Start update_ice_slow_thermo", FIA, sG, US)
! call IOF_chksum("Start update_ice_slow_thermo", Ice%sCS%IOF, sG, US)
endif
! Store some diagnostic fluxes...
!$OMP parallel do default(none) shared(sG, FIA)
do j=sG%jsc,sG%jec ; do i=sG%isc,sG%iec
FIA%calving_preberg(i,j) = FIA%calving(i,j)
FIA%calving_hflx_preberg(i,j) = FIA%calving_hflx(i,j)
enddo ; enddo
if (Ice%sCS%redo_fast_update) then
call redo_update_ice_model_fast(sIST, Ice%sCS%sOSS, Ice%sCS%Rad, FIA, Ice%sCS%TSF, &
Ice%sCS%optics_CSp, Ice%sCS%Time_step_slow, Ice%sCS%fast_thermo_CSp, sG, US, sIG)
call find_excess_fluxes(FIA, Ice%sCS%TSF, Ice%sCS%XSF, sIST%part_size, sG, US, sIG)
endif
call convert_frost_to_snow(FIA, sG, US, sIG)
if (Ice%sCS%do_icebergs) then
if (Ice%sCS%berg_windstress_bug) then
! This code is only required to reproduce an old bug.
i_off = LBOUND(Ice%flux_t,1) - sG%isc
j_off = LBOUND(Ice%flux_t,2) - sG%jsc
!$OMP parallel do default(none) shared(Ice,sG,US,i_off,j_off) private(i2,j2)
do j=sG%jsc,sG%jec ; do i=sG%isc,sG%iec
i2 = i+i_off ; j2 = j+j_off
Ice%sCS%IOF%flux_u_ocn(i,j) = US%kg_m2s_to_RZ_T*US%m_s_to_L_T*Ice%flux_u(i2,j2)
Ice%sCS%IOF%flux_v_ocn(i,j) = US%kg_m2s_to_RZ_T*US%m_s_to_L_T*Ice%flux_v(i2,j2)
enddo ; enddo
endif
call cpu_clock_end(ice_clock_slow) ; call cpu_clock_end(iceClock)
call update_icebergs(sIST, Ice%sCS%OSS, Ice%sCS%IOF, FIA, Ice%icebergs, US%T_to_s*dt_slow, &
sG, US, sIG, Ice%sCS%dyn_trans_CSp)
call cpu_clock_begin(iceClock) ; call cpu_clock_begin(ice_clock_slow)
if (Ice%sCS%debug) then
call FIA_chksum("After update_icebergs", FIA, sG, US)
endif
endif
if (Ice%sCS%debug) then
call Ice_public_type_chksum("Before slow_thermodynamics", Ice, check_slow=.true.)
call FIA_chksum("Before slow_thermodynamics", FIA, sG, US)
call IST_chksum("Before slow_thermodynamics", sIST, sG, US, sIG)
call OSS_chksum("Before slow_thermodynamics", Ice%sCS%OSS, sG, US)
if (associated(Ice%sCS%XSF)) &
call TSF_chksum("Before slow_thermodynamics XSF", Ice%sCS%XSF, sG, US)
! call IOF_chksum("Before slow_thermodynamics", Ice%sCS%IOF, sG, US)
endif
call slow_thermodynamics(sIST, dt_slow, Ice%sCS%slow_thermo_CSp, Ice%sCS%OSS, FIA, &
Ice%sCS%XSF, Ice%sCS%IOF, sG, US, sIG)
if (Ice%sCS%debug) then
call Ice_public_type_chksum("Before set_ocean_top_fluxes", Ice, check_slow=.true.)
call IOF_chksum("Before set_ocean_top_fluxes", Ice%sCS%IOF, sG, US, thermo_fluxes=.true.)
call IST_chksum("Before set_ocean_top_fluxes", sIST, sG, US, sIG)
endif
! Set up the thermodynamic fluxes in the externally visible structure Ice.
call set_ocean_top_fluxes(Ice, sIST, Ice%sCS%IOF, FIA, Ice%sCS%OSS, sG, US, sIG, Ice%sCS)
call cpu_clock_end(ice_clock_slow) ; call cpu_clock_end(iceClock)
end subroutine update_ice_slow_thermo
!-----------------------------------------------------------------------
!> Update the sea-ice state due to dynamics and ice transport.
subroutine update_ice_dynamics_trans(Ice, time_step, start_cycle, end_cycle, cycle_length)
type(ice_data_type), intent(inout) :: Ice !< The publicly visible ice data type.
type(time_type), optional, intent(in) :: time_step !< The amount of time to cover in this update.
logical, optional, intent(in) :: start_cycle !< This indicates whether this call is to be
!! treated as the first call to update_ice_dynamics_trans
!! in a time-stepping cycle; missing is like true.
logical, optional, intent(in) :: end_cycle !< This indicates whether this call is to be
!! treated as the last call to update_ice_dynamics_trans
!! in a time-stepping cycle; missing is like true.
real, optional, intent(in) :: cycle_length !< The duration of a coupled time stepping cycle [s].
! These pointers are used to simplify the code below.
type(ice_grid_type), pointer :: sIG => NULL()
type(SIS_hor_grid_type), pointer :: sG => NULL()
type(ice_state_type), pointer :: sIST => NULL()
type(fast_ice_avg_type), pointer :: FIA => NULL()
type(unit_scale_type), pointer :: US => NULL()
real :: dt_slow ! The time step over which to advance the model [T ~> s].
logical :: do_multi_trans, cycle_start
if (.not.associated(Ice%sCS)) call SIS_error(FATAL, &
"The pointer to Ice%sCS must be associated in update_ice_dynamics_trans.")
sIST => Ice%sCS%IST ; sIG => Ice%sCS%IG ; sG => Ice%sCS%G ; FIA => Ice%sCS%FIA ; US => Ice%sCS%US
dt_slow = US%s_to_T*time_type_to_real(Ice%sCS%Time_step_slow)
if (present(time_step)) dt_slow = US%s_to_T*time_type_to_real(time_step)
cycle_start = .true. ; if (present(start_cycle)) cycle_start = start_cycle
call cpu_clock_begin(iceClock) ; call cpu_clock_begin(ice_clock_slow)
! Do halo updates on the forcing fields, as necessary. This must occur before
! the call to SIS_dynamics_trans, because update_icebergs does its own halo
! updates, and slow_thermodynamics only works on the computational domain.
if (cycle_start) then
call pass_vector(FIA%WindStr_x, FIA%WindStr_y, sG%Domain, stagger=AGRID, complete=.false.)
call pass_vector(FIA%WindStr_ocn_x, FIA%WindStr_ocn_y, sG%Domain, stagger=AGRID)
call pass_var(FIA%ice_cover, sG%Domain, complete=.false.)
call pass_var(FIA%ice_free, sG%Domain, complete=.true.)
endif
if (sIST%valid_IST) then
call pass_var(sIST%part_size, sG%Domain)
call pass_var(sIST%mH_ice, sG%Domain, complete=.false.)
call pass_var(sIST%mH_pond, sG%Domain, complete=.false.)
call pass_var(sIST%mH_snow, sG%Domain, complete=.true.)
endif
if (Ice%sCS%debug) then
call Ice_public_type_chksum("Before SIS_dynamics_trans", Ice, check_slow=.true.)
endif
do_multi_trans = (present(start_cycle) .or. present(end_cycle) .or. present(cycle_length))
if (Ice%sCS%specified_ice) then ! There is no ice dynamics or transport.
call specified_ice_dynamics(sIST, Ice%sCS%OSS, FIA, Ice%sCS%IOF, dt_slow, &
Ice%sCS%specified_ice_CSp, sG, US, sIG)
elseif (do_multi_trans) then
call SIS_multi_dyn_trans(sIST, Ice%sCS%OSS, FIA, Ice%sCS%IOF, dt_slow, Ice%sCS%dyn_trans_CSp, &
Ice%icebergs, sG, US, sIG, Ice%sCS%SIS_tracer_flow_CSp, &
Ice%OBC, start_cycle, end_cycle, cycle_length)
elseif (Ice%sCS%slab_ice) then ! Use a very old slab ice model.
call slab_ice_dyn_trans(sIST, Ice%sCS%OSS, FIA, Ice%sCS%IOF, dt_slow, Ice%sCS%dyn_trans_CSp, &
sG, US, sIG, Ice%sCS%SIS_tracer_flow_CSp, Ice%OBC)
else ! This is the typical branch used by SIS2.
call SIS_dynamics_trans(sIST, Ice%sCS%OSS, FIA, Ice%sCS%IOF, dt_slow, Ice%sCS%dyn_trans_CSp, &
Ice%icebergs, sG, US, sIG, Ice%sCS%SIS_tracer_flow_CSp, Ice%OBC)
endif
! Set up the stresses and surface pressure in the externally visible structure Ice.
if (sIST%valid_IST) call ice_mass_from_IST(sIST, Ice%sCS%IOF, sG, sIG)
if (Ice%sCS%debug) then
call IOF_chksum("Before set_ocean_top_dyn_fluxes", Ice%sCS%IOF, sG, US, mech_fluxes=.true.)
endif
call set_ocean_top_dyn_fluxes(Ice, Ice%sCS%IOF, FIA, sG, US, Ice%sCS)
if (Ice%sCS%debug) then
call Ice_public_type_chksum("End update_ice_dynamics_trans", Ice, check_slow=.true.)
endif
!### THIS NO LONGER WORKS ON SLOW ICE PES.
! if (Ice%sCS%bounds_check) then
! call Ice_public_type_bounds_check(Ice, sG, "End update_ice_slow")
! endif
call cpu_clock_end(ice_clock_slow) ; call cpu_clock_end(iceClock)
end subroutine update_ice_dynamics_trans
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
!> ice_model_fast_cleanup performs the final steps in the fast ice update cycle
!! and prepares data to drive the slow ice updates. This includes finding the
!! averaged fluxes and unpacking the land to ice forcing.
subroutine ice_model_fast_cleanup(Ice)
type(ice_data_type), intent(inout) :: Ice !< The publicly visible ice data type.
if (.not.associated(Ice%fCS)) call SIS_error(FATAL, &
"The pointer to Ice%fCS must be associated in ice_model_fast_cleanup.")
! average fluxes from update_ice_model_fast
call avg_top_quantities(Ice%fCS%FIA, Ice%fCS%Rad, Ice%fCS%IST, Ice%fCS%G, Ice%fCS%US, Ice%fCS%IG)
call total_top_quantities(Ice%fCS%FIA, Ice%fCS%TSF, Ice%fCS%IST%part_size, Ice%fCS%G, Ice%fCS%US, Ice%fCS%IG)
if (allocated(Ice%fCS%IST%t_surf)) &
Ice%fCS%IST%t_surf(:,:,1:) = Ice%fCS%Rad%T_skin(:,:,:) + Ice%fCS%IST%T_0degC
call infill_array(Ice%fCS%IST, Ice%fCS%sOSS%T_fr_ocn, Ice%fCS%Rad%T_skin, Ice%fCS%G, Ice%fCS%IG)
end subroutine ice_model_fast_cleanup
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
!> unpack_land_ice_bdry converts the information in a publicly visible
!! land_ice_boundary_type into an internally visible fast_ice_avg_type variable.
subroutine unpack_land_ice_boundary(Ice, LIB)
type(ice_data_type), intent(inout) :: Ice !< The publicly visible ice data type.
type(land_ice_boundary_type), intent(in) :: LIB !< The land ice boundary type that is being unpacked.
type(fast_ice_avg_type), pointer :: FIA => NULL()
type(SIS_hor_grid_type), pointer :: G => NULL()
type(unit_scale_type), pointer :: US => NULL()
integer :: i, j, k, m, n, i2, j2, k2, isc, iec, jsc, jec, i_off, j_off
if (.not.associated(Ice%fCS)) call SIS_error(FATAL, &
"The pointer to Ice%fCS must be associated in unpack_land_ice_boundary.")
if (.not.associated(Ice%fCS%FIA)) call SIS_error(FATAL, &
"The pointer to Ice%fCS%FIA must be associated in unpack_land_ice_boundary.")
if (.not.associated(Ice%fCS%G)) call SIS_error(FATAL, &
"The pointer to Ice%fCS%G must be associated in unpack_land_ice_boundary.")
FIA => Ice%fCS%FIA ; G => Ice%fCS%G
US => Ice%fCS%US
isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec
! Store liquid runoff and other fluxes from the land to the ice or ocean.
i_off = LBOUND(LIB%runoff,1) - G%isc ; j_off = LBOUND(LIB%runoff,2) - G%jsc
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,FIA,LIB,i_off,j_off,G,US) &
!$OMP private(i2,j2)
do j=jsc,jec ; do i=isc,iec ; if (G%mask2dT(i,j) > 0.0) then
i2 = i+i_off ; j2 = j+j_off
FIA%runoff(i,j) = US%kg_m2s_to_RZ_T*LIB%runoff(i2,j2)
FIA%calving(i,j) = US%kg_m2s_to_RZ_T*LIB%calving(i2,j2)
FIA%runoff_hflx(i,j) = US%W_m2_to_QRZ_T*LIB%runoff_hflx(i2,j2)
FIA%calving_hflx(i,j) = US%W_m2_to_QRZ_T*LIB%calving_hflx(i2,j2)
else
! This is a land point from the perspective of the sea-ice.
! At some point it might make sense to check for non-zero fluxes, which
! might indicate regridding errors. However, bad-data values are also
! non-zero and should not be flagged.
FIA%runoff(i,j) = 0.0 ; FIA%calving(i,j) = 0.0
FIA%runoff_hflx(i,j) = 0.0 ; FIA%calving_hflx(i,j) = 0.0
endif ; enddo ; enddo
if (Ice%fCS%debug) then
call FIA_chksum("End of unpack_land_ice_boundary", FIA, G, Ice%fCS%US)
endif
end subroutine unpack_land_ice_boundary
!> unpack_ocean_ice_boundary_calved_shelf_bergs converts the calving information in a publicly visible
!! ocean_ice_boundary_type into an internally visible fast_ice_avg_type variable.
subroutine unpack_ocean_ice_boundary_calved_shelf_bergs(Ice, OIB)
type(ice_data_type), intent(inout) :: Ice !< The publicly visible ice data type.
type(ocean_ice_boundary_type), intent(in) :: OIB !< The ocean ice boundary type that is being unpacked.
type(fast_ice_avg_type), pointer :: FIA => NULL()
type(SIS_hor_grid_type), pointer :: G => NULL()
type(unit_scale_type), pointer :: US => NULL()
integer :: i, j, k, m, n, i2, j2, k2, isc, iec, jsc, jec, i_off, j_off
if (.not.associated(Ice%fCS)) call SIS_error(FATAL, &
"The pointer to Ice%fCS must be associated in unpack_ocean_ice_boundary_calved_shelf_bergs.")
if (.not.associated(Ice%fCS%FIA)) call SIS_error(FATAL, &
"The pointer to Ice%fCS%FIA must be associated in unpack_ocean_ice_boundary_calved_shelf_berg.")
if (.not.associated(Ice%fCS%G)) call SIS_error(FATAL, &
"The pointer to Ice%fCS%G must be associated in unpack_ocean_ice_boundary_calved_shelf_berg.")
FIA => Ice%fCS%FIA ; G => Ice%fCS%G
US => Ice%fCS%US
isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec
! Store calving flux from ice shelves to the sea ice or ocean.
i_off = LBOUND(OIB%calving,1) - G%isc ; j_off = LBOUND(OIB%calving,2) - G%jsc
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,FIA,OIB,i_off,j_off,G,US) &
!$OMP private(i2,j2)
do j=jsc,jec ; do i=isc,iec ; if (G%mask2dT(i,j) > 0.0) then
i2 = i+i_off ; j2 = j+j_off
if (OIB%calving(i2,j2)>0.0) then
if (FIA%calving(i,j)>0.0) call SIS_error(FATAL,"Overlap in calving from snow discharge and ice shelf!")
FIA%calving(i,j) = US%kg_m2s_to_RZ_T*OIB%calving(i2,j2)
FIA%calving_hflx(i,j) = US%W_m2_to_QRZ_T*OIB%calving_hflx(i2,j2)
endif
endif ; enddo ; enddo
if (Ice%fCS%debug) then
call FIA_chksum("End of unpack_ocean_ice_boundary_calved_shelf_berg", FIA, G, Ice%fCS%US)
endif
end subroutine unpack_ocean_ice_boundary_calved_shelf_bergs
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
!> This subroutine copies information (mostly fluxes and the updated temperatures)
!! from the fast part of the sea-ice to the slow part of the sea ice.
subroutine exchange_fast_to_slow_ice(Ice)
type(ice_data_type), &
intent(inout) :: Ice !< The publicly visible ice data type whose fast
!! part is to be exchanged with the slow part.
type(fast_ice_avg_type), pointer :: FIA_null => NULL()
type(ice_state_type), pointer :: IST_null => NULL()
type(ice_rad_type), pointer :: Rad_null => NULL()
type(total_sfc_flux_type), pointer :: TSF_null => NULL()
integer :: isc, iec, jsc, jec, isd, ied, jsd, jed
logical :: redo_fast_update
redo_fast_update = .false.
if (associated(Ice%fCS)) redo_fast_update = Ice%fCS%redo_fast_update
if (associated(Ice%sCS)) redo_fast_update = Ice%sCS%redo_fast_update
if (associated(Ice%fCS)) then
isc = Ice%fCS%G%isc ; iec = Ice%fCS%G%iec ; jsc = Ice%fCS%G%jsc ; jec = Ice%fCS%G%jec
isd = Ice%fCS%G%isd ; ied = Ice%fCS%G%ied ; jsd = Ice%fCS%G%jsd ; jed = Ice%fCS%G%jed
! Propagate the coupler_type info to Ice%fCS%FIA%tr_flux and allocate its arrays.
call coupler_type_spawn(Ice%ocean_fluxes, Ice%fCS%FIA%tr_flux, &
(/isd, isc, iec, ied/), (/jsd, jsc, jec, jed/), &
(/0, Ice%fCS%IG%CatIce/), as_needed=.true.)
if (redo_fast_update) &
! Propagate the coupler_type info to Ice%fCS%TSF%tr_flux and allocate its arrays.
call coupler_type_spawn(Ice%ocean_fluxes, Ice%fCS%TSF%tr_flux, &
(/isd, isc, iec, ied/), (/jsd, jsc, jec, jed/), as_needed=.true.)
endif
if (associated(Ice%sCS)) then
isc = Ice%sCS%G%isc ; iec = Ice%sCS%G%iec ; jsc = Ice%sCS%G%jsc ; jec = Ice%sCS%G%jec
isd = Ice%sCS%G%isd ; ied = Ice%sCS%G%ied ; jsd = Ice%sCS%G%jsd ; jed = Ice%sCS%G%jed
! Propagate the coupler_type info to Ice%sCS%FIA%tr_flux and allocate its arrays.
call coupler_type_spawn(Ice%ocean_fluxes, Ice%sCS%FIA%tr_flux, &
(/isd, isc, iec, ied/), (/jsd, jsc, jec, jed/), &
(/0, Ice%sCS%IG%CatIce/), as_needed=.true.)
if (redo_fast_update) &
! Propagate the coupler_type info to Ice%sCS%TSF%tr_flux and allocate its arrays.
call coupler_type_spawn(Ice%ocean_fluxes, Ice%sCS%TSF%tr_flux, &
(/isd, isc, iec, ied/), (/jsd, jsc, jec, jed/), as_needed=.true.)
endif
if (Ice%xtype == DIRECT) then
if (.not.associated(Ice%fCS) .or. .not.associated(Ice%sCS)) call SIS_error(FATAL, &
"With xtype=DIRECT, both the pointer to Ice%sCS and the pointer to Ice%fCS must be "//&
"associated (although perhaps not with each other) in exchange_fast_to_slow_ice.")
if (.not.associated(Ice%fCS%FIA, Ice%sCS%FIA)) then
call copy_FIA_to_FIA(Ice%fCS%FIA, Ice%sCS%FIA, Ice%fCS%G%HI, Ice%sCS%G%HI, Ice%sCS%IG)
endif
if (redo_fast_update) then
if (.not.associated(Ice%fCS%TSF, Ice%sCS%TSF)) &
call copy_TSF_to_TSF(Ice%fCS%TSF, Ice%sCS%TSF, Ice%fCS%G%HI, Ice%sCS%G%HI)
if (.not.associated(Ice%fCS%Rad, Ice%sCS%Rad)) &
call copy_Rad_to_Rad(Ice%fCS%Rad, Ice%sCS%Rad, Ice%fCS%G%HI, Ice%sCS%G%HI, Ice%fCS%IG)
else
if (.not.associated(Ice%fCS%IST, Ice%sCS%IST)) &
call copy_IST_to_IST(Ice%fCS%IST, Ice%sCS%IST, Ice%fCS%G%HI, Ice%sCS%G%HI, Ice%fCS%IG)
endif
elseif (Ice%xtype == REDIST) then
if (.not.associated(Ice%fCS) .and. .not.associated(Ice%sCS)) call SIS_error(FATAL, &
"Either the pointer to Ice%sCS or the pointer to Ice%fCS must be "//&
"associated in exchange_fast_to_slow_ice.")
if (associated(Ice%fCS) .and. associated(Ice%sCS)) then
if (.not.associated(Ice%fCS%FIA, Ice%sCS%FIA)) &
call redistribute_FIA_to_FIA(Ice%fCS%FIA, Ice%sCS%FIA, Ice%fast_domain, &
Ice%slow_domain, Ice%sCS%G, Ice%sCS%IG)
if (redo_fast_update) then
call redistribute_TSF_to_TSF(Ice%fCS%TSF, Ice%sCS%TSF, Ice%fast_domain, &
Ice%slow_domain, Ice%sCS%G%HI)
call redistribute_Rad_to_Rad(Ice%fCS%Rad, Ice%sCS%Rad, Ice%fast_domain, Ice%slow_domain)
else
if (.not.associated(Ice%fCS%IST, Ice%sCS%IST)) &
call redistribute_IST_to_IST(Ice%fCS%IST, Ice%sCS%IST, Ice%fast_domain, Ice%slow_domain)
endif
elseif (associated(Ice%fCS)) then
call redistribute_FIA_to_FIA(Ice%fCS%FIA, FIA_null, Ice%fast_domain, Ice%slow_domain)
if (redo_fast_update) then
call redistribute_TSF_to_TSF(Ice%fCS%TSF, TSF_null, Ice%fast_domain, Ice%slow_domain)
call redistribute_Rad_to_Rad(Ice%fCS%Rad, Rad_null, Ice%fast_domain, Ice%slow_domain)
else
call redistribute_IST_to_IST(Ice%fCS%IST, IST_null, Ice%fast_domain, Ice%slow_domain)
endif
elseif (associated(Ice%sCS)) then
call redistribute_FIA_to_FIA(FIA_null, Ice%sCS%FIA, Ice%fast_domain, &
Ice%slow_domain, Ice%sCS%G, Ice%sCS%IG)
if (redo_fast_update) then
call redistribute_TSF_to_TSF(TSF_null, Ice%sCS%TSF, Ice%fast_domain, &
Ice%slow_domain, Ice%sCS%G%HI)
call redistribute_Rad_to_Rad(Rad_null, Ice%sCS%Rad, Ice%fast_domain, Ice%slow_domain)
else
call redistribute_IST_to_IST(IST_null, Ice%sCS%IST, Ice%fast_domain, Ice%slow_domain)
endif
else
call SIS_error(FATAL, "Either the pointer to Ice%sCS or the pointer to "//&
"Ice%fCS must be associated in exchange_fast_to_slow_ice.")
endif
else
call SIS_error(FATAL, "exchange_fast_to_slow_ice called with an unrecognized Ice%xtype value.")
endif
end subroutine exchange_fast_to_slow_ice
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
!> set_ocean_top_fluxes translates ice-bottom fluxes of heat, mass, salt, and
!! tracers from the ice model's internal state to the public ice data type
!! for use by the ocean model.
subroutine set_ocean_top_fluxes(Ice, IST, IOF, FIA, OSS, G, US, IG, sCS)
type(ice_data_type), intent(inout) :: Ice !< The publicly visible ice data type.
type(ice_state_type), intent(inout) :: IST !< A type describing the state of the sea ice
type(ice_ocean_flux_type), intent(in) :: IOF !< A structure containing fluxes from the ice to
!! the ocean that are calculated by the ice model.
type(fast_ice_avg_type), intent(in) :: FIA !< A type containing averages of fields
!! (mostly fluxes) over the fast updates
type(ocean_sfc_state_type), intent(in) :: OSS !< A structure containing the arrays that describe
!! the ocean's surface state for the ice model.
type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type
type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors
type(ice_grid_type), intent(in) :: IG !< The sea-ice specific grid type
type(SIS_slow_CS), intent(in) :: sCS !< The slow ice control structure
real :: I_count
integer :: i, j, k, isc, iec, jsc, jec, m, n
integer :: i2, j2, i_off, j_off, ncat, NkIce
isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec
ncat = IG%CatIce ; NkIce = IG%NkIce
if (sCS%debug) then
call Ice_public_type_chksum("Start set_ocean_top_fluxes", Ice, check_slow=.true.)
call IOF_chksum("Start set_ocean_top_fluxes", IOF, G, sCS%US, thermo_fluxes=.true.)
call FIA_chksum("Start set_ocean_top_fluxes", FIA, G, US)
endif
! It is possible that the ice mass and surface pressure will be needed after
! the thermodynamic step, in which case this should be uncommented.
! ! Sum the concentration weighted mass.
! Ice%mi(:,:) = 0.0
! i_off = LBOUND(Ice%mi,1) - G%isc ; j_off = LBOUND(Ice%mi,2) - G%jsc
! !$OMP parallel do default(shared) private(i2,j2)
! do j=jsc,jec ; do k=1,ncat ; do i=isc,iec
! i2 = i+i_off ; j2 = j+j_off! Use these to correct for indexing differences.
! Ice%mi(i2,j2) = Ice%mi(i2,j2) + IST%part_size(i,j,k) * &
! (G%US%RZ_to_kg_m2 * ((IST%mH_snow(i,j,k) + IST%mH_pond(i,j,k)) + IST%mH_ice(i,j,k)))
! enddo ; enddo ; enddo
! This block of code is probably unnecessary.
Ice%flux_t(:,:) = 0.0 ; Ice%flux_q(:,:) = 0.0
Ice%flux_sw_nir_dir(:,:) = 0.0 ; Ice%flux_sw_nir_dif(:,:) = 0.0
Ice%flux_sw_vis_dir(:,:) = 0.0 ; Ice%flux_sw_vis_dif(:,:) = 0.0
Ice%flux_lw(:,:) = 0.0 ; Ice%flux_lh(:,:) = 0.0
Ice%fprec(:,:) = 0.0 ; Ice%lprec(:,:) = 0.0
call coupler_type_rescale_data(Ice%ocean_fluxes, 0.0)
i_off = LBOUND(Ice%flux_t,1) - G%isc ; j_off = LBOUND(Ice%flux_t,2) - G%jsc
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,Ice,IST,IOF,FIA,i_off,j_off,G,US,OSS) &
!$OMP private(i2,j2)
do j=jsc,jec ; do i=isc,iec
i2 = i+i_off ; j2 = j+j_off! Use these to correct for indexing differences.
Ice%flux_t(i2,j2) = US%QRZ_T_to_W_m2*IOF%flux_sh_ocn_top(i,j)
Ice%flux_q(i2,j2) = US%RZ_T_to_kg_m2s*IOF%evap_ocn_top(i,j)
Ice%flux_sw_vis_dir(i2,j2) = US%QRZ_T_to_W_m2*IOF%flux_sw_ocn(i,j,VIS_DIR)
Ice%flux_sw_vis_dif(i2,j2) = US%QRZ_T_to_W_m2*IOF%flux_sw_ocn(i,j,VIS_DIF)
Ice%flux_sw_nir_dir(i2,j2) = US%QRZ_T_to_W_m2*IOF%flux_sw_ocn(i,j,NIR_DIR)
Ice%flux_sw_nir_dif(i2,j2) = US%QRZ_T_to_W_m2*IOF%flux_sw_ocn(i,j,NIR_DIF)
Ice%flux_lw(i2,j2) = US%QRZ_T_to_W_m2*IOF%flux_lw_ocn_top(i,j)
Ice%flux_lh(i2,j2) = US%QRZ_T_to_W_m2*IOF%flux_lh_ocn_top(i,j)
Ice%fprec(i2,j2) = US%RZ_T_to_kg_m2s*IOF%fprec_ocn_top(i,j)
Ice%lprec(i2,j2) = US%RZ_T_to_kg_m2s*IOF%lprec_ocn_top(i,j)
Ice%runoff(i2,j2) = US%RZ_T_to_kg_m2s*FIA%runoff(i,j)
Ice%calving(i2,j2) = US%RZ_T_to_kg_m2s*FIA%calving(i,j)
Ice%runoff_hflx(i2,j2) = US%QRZ_T_to_W_m2*FIA%runoff_hflx(i,j)
Ice%calving_hflx(i2,j2) = US%QRZ_T_to_W_m2*FIA%calving_hflx(i,j)
Ice%flux_salt(i2,j2) = US%S_to_ppt*US%RZ_T_to_kg_m2s*IOF%flux_salt(i,j)
Ice%SST_C(i2,j2) = US%C_to_degC*OSS%SST_C(i,j)
! It is possible that the ice mass and surface pressure will be needed after
! the thermodynamic step, in which case this should be uncommented.
! if (IOF%slp2ocean) then
! Ice%p_surf(i2,j2) = US%RZ_T_to_kg_m2s*US%L_T_to_m_s*FIA%p_atm_surf(i,j) - 1e5 ! SLP - 1 std. atmosphere [Pa].
! else
! Ice%p_surf(i2,j2) = 0.0
! endif
! Ice%p_surf(i2,j2) = Ice%p_surf(i2,j2) + US%L_T_to_m_s**2*US%m_to_Z*G%g_Earth*Ice%mi(i2,j2)
enddo ; enddo
if (allocated(IOF%melt_nudge)) then
do j=jsc,jec ; do i=isc,iec
i2 = i+i_off ; j2 = j+j_off! Use these to correct for indexing differences.
Ice%lprec(i2,j2) = Ice%lprec(i2,j2) + US%RZ_T_to_kg_m2s*IOF%melt_nudge(i,j)
enddo ; enddo
endif
! This copy may need to be skipped in the first step of a cold-start run with lagged ice
! coupling, but otherwise if it is skipped may indicate a problem that should be trapped.
if (coupler_type_initialized(IOF%tr_flux_ocn_top)) &
call coupler_type_copy_data(IOF%tr_flux_ocn_top, Ice%ocean_fluxes)
if (sCS%debug) then
call Ice_public_type_chksum("End set_ocean_top_fluxes", Ice, check_slow=.true.)
endif
end subroutine set_ocean_top_fluxes
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
!> ice_mass_from_IST stores the total ice mass determined from IST in the IOF type.
subroutine ice_mass_from_IST(IST, IOF, G, IG)
type(ice_state_type), intent(inout) :: IST !< A type describing the state of the sea ice
type(ice_ocean_flux_type), intent(inout) :: IOF !< A structure containing fluxes from the ice to
!! the ocean that are calculated by the ice model.
type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type
type(ice_grid_type), intent(in) :: IG !< The sea-ice specific grid type
integer :: i, j, k, isc, iec, jsc, jec, ncat
isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; ncat = IG%CatIce
! Sum the concentration weighted mass.
IOF%mass_ice_sn_p(:,:) = 0.0
!$OMP parallel do default(shared)
do j=jsc,jec ; do k=1,ncat ; do i=isc,iec
IOF%mass_ice_sn_p(i,j) = IOF%mass_ice_sn_p(i,j) + IST%part_size(i,j,k) * &
((IST%mH_snow(i,j,k) + IST%mH_pond(i,j,k)) + IST%mH_ice(i,j,k))
enddo ; enddo ; enddo
end subroutine ice_mass_from_IST
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
!> set_ocean_top_dyn_fluxes translates ice-bottom stresses and mass from the ice
!! model's ice-ocean flux type and the fast-ice average type to the public
!! ice data type for use by the ocean model.
subroutine set_ocean_top_dyn_fluxes(Ice, IOF, FIA, G, US, sCS)
type(ice_data_type), intent(inout) :: Ice !< The publicly visible ice data type.
type(ice_ocean_flux_type), intent(in) :: IOF !< A structure containing fluxes from the ice to
!! the ocean that are calculated by the ice model.
type(fast_ice_avg_type), intent(in) :: FIA !< A type containing averages of fields
!! (mostly fluxes) over the fast updates
type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type
type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors
type(SIS_slow_CS), intent(in) :: sCS !< The slow ice control structure
real :: I_count
integer :: i, j, k, isc, iec, jsc, jec
integer :: i2, j2, i_off, j_off, ind
isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec
if (sCS%debug) then
call Ice_public_type_chksum("Start set_ocean_top_dyn_fluxes", Ice, check_slow=.true.)
call IOF_chksum("Start set_ocean_top_dyn_fluxes", IOF, G, US, mech_fluxes=.true.)
endif
! Sum the concentration weighted mass.
Ice%mi(:,:) = 0.0
i_off = LBOUND(Ice%mi,1) - G%isc ; j_off = LBOUND(Ice%mi,2) - G%jsc
!$OMP parallel do default(shared) private(i2,j2)
do j=jsc,jec ; do i=isc,iec
i2 = i+i_off ; j2 = j+j_off! Use these to correct for indexing differences.
Ice%mi(i2,j2) = Ice%mi(i2,j2) + US%RZ_to_kg_m2*IOF%mass_ice_sn_p(i,j)
enddo ; enddo
if (sCS%do_icebergs .and. associated(IOF%mass_berg)) then
! Note that the IOF berg fields and Ice fields are only allocated on the
! computational domains, although they may use different indexing conventions.
Ice%mi(:,:) = Ice%mi(:,:) + IOF%mass_berg(:,:)
if (sCS%pass_iceberg_area_to_ocean) then
Ice%mass_berg(:,:) = IOF%mass_berg(:,:)
if (associated(IOF%ustar_berg)) Ice%ustar_berg(:,:) = IOF%ustar_berg(:,:)
if (associated(IOF%area_berg)) Ice%area_berg(:,:) = IOF%area_berg(:,:)
endif
endif
i_off = LBOUND(Ice%flux_t,1) - G%isc ; j_off = LBOUND(Ice%flux_t,2) - G%jsc
!$OMP parallel do default(shared) private(i2,j2)
do j=jsc,jec ; do i=isc,iec
i2 = i+i_off ; j2 = j+j_off! Use these to correct for indexing differences.
Ice%flux_u(i2,j2) = US%RZ_T_to_kg_m2s*US%L_T_to_m_s*IOF%flux_u_ocn(i,j)
Ice%flux_v(i2,j2) = US%RZ_T_to_kg_m2s*US%L_T_to_m_s*IOF%flux_v_ocn(i,j)
if (IOF%slp2ocean) then
Ice%p_surf(i2,j2) = US%RZ_T_to_kg_m2s*US%L_T_to_m_s*FIA%p_atm_surf(i,j) - 1e5 ! SLP - 1 std. atmosphere [Pa].
else
Ice%p_surf(i2,j2) = 0.0
endif
Ice%p_surf(i2,j2) = Ice%p_surf(i2,j2) + US%L_T_to_m_s**2*US%m_to_Z*G%g_Earth*Ice%mi(i2,j2)
enddo ; enddo
if (associated(Ice%stress_mag) .and. allocated(IOF%stress_mag)) then
i_off = LBOUND(Ice%stress_mag,1) - G%isc ; j_off = LBOUND(Ice%stress_mag,2) - G%jsc
!$OMP parallel do default(shared) private(i2,j2)
do j=jsc,jec ; do i=isc,iec ; i2 = i+i_off ; j2 = j+j_off
Ice%stress_mag(i2,j2) = US%RZ_T_to_kg_m2s*US%L_T_to_m_s*IOF%stress_mag(i,j)
enddo ; enddo
endif
if (sCS%debug) then
call Ice_public_type_chksum("End set_ocean_top_dyn_fluxes", Ice, check_slow=.true.)
endif
end subroutine set_ocean_top_dyn_fluxes
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
!> This subroutine copies information from the slow part of the sea-ice to the
!! fast part of the sea ice.
subroutine exchange_slow_to_fast_ice(Ice)
type(ice_data_type), &
intent(inout) :: Ice !< The publicly visible ice data type whose slow
!! part is to be exchanged with the fast part.
type(simple_OSS_type), pointer :: sOSS_null => NULL()
type(ice_state_type), pointer :: IST_null => NULL()
integer :: isc, iec, jsc, jec, isd, ied, jsd, jed
call cpu_clock_begin(iceClock) ; call cpu_clock_begin(ice_clock_exchange)
if (associated(Ice%fCS)) then
isc = Ice%fCS%G%isc ; iec = Ice%fCS%G%iec ; jsc = Ice%fCS%G%jsc ; jec = Ice%fCS%G%jec
isd = Ice%fCS%G%isd ; ied = Ice%fCS%G%ied ; jsd = Ice%fCS%G%jsd ; jed = Ice%fCS%G%jed
! Propagate the coupler_type info to Ice%fCS%sOSS%tr_fields and allocate its arrays.
call coupler_type_spawn(Ice%ocean_fields, Ice%fCS%sOSS%tr_fields, &
(/isd, isc, iec, ied/), (/jsd, jsc, jec, jed/), as_needed=.true. )
endif
if (Ice%xtype == DIRECT) then
if (.not.associated(Ice%fCS) .or. .not.associated(Ice%sCS)) call SIS_error(FATAL, &
"With xtype=DIRECT, both the pointer to Ice%sCS and the pointer to Ice%fCS must be "//&
"associated (although perhaps not with each other) in exchange_slow_to_fast_ice.")
if (.not.associated(Ice%fCS%sOSS, Ice%sCS%sOSS)) then
call copy_sOSS_to_sOSS(Ice%sCS%sOSS, Ice%fCS%sOSS, Ice%sCS%G%HI, Ice%fCS%G%HI)
endif
if (.not.associated(Ice%fCS%IST, Ice%sCS%IST)) then
call copy_IST_to_IST(Ice%sCS%IST, Ice%fCS%IST, Ice%sCS%G%HI, Ice%fCS%G%HI, &
Ice%sCS%IG)
endif
elseif (Ice%xtype == REDIST) then
if (.not.associated(Ice%fCS) .and. .not.associated(Ice%sCS)) call SIS_error(FATAL, &
"Either the pointer to Ice%sCS or the pointer to Ice%fCS must be "//&
"associated in exchange_slow_to_fast_ice.")
if (associated(Ice%fCS) .and. associated(Ice%sCS)) then
if (.not.associated(Ice%fCS%sOSS, Ice%sCS%sOSS)) &
call redistribute_sOSS_to_sOSS(Ice%sCS%sOSS, Ice%fCS%sOSS, Ice%slow_domain, &
Ice%fast_domain, Ice%fCS%G%HI)
if (.not.associated(Ice%fCS%IST, Ice%sCS%IST)) &
call redistribute_IST_to_IST(Ice%sCS%IST, Ice%fCS%IST, Ice%slow_domain, &
Ice%fast_domain)
elseif (associated(Ice%fCS)) then
call redistribute_sOSS_to_sOSS(sOSS_null, Ice%fCS%sOSS, Ice%slow_domain, &
Ice%fast_domain, HI_out=Ice%fCS%G%HI)
call redistribute_IST_to_IST(IST_null, Ice%fCS%IST, Ice%slow_domain, &
Ice%fast_domain)
elseif (associated(Ice%sCS)) then
call redistribute_sOSS_to_sOSS(Ice%sCS%sOSS, sOSS_null, Ice%slow_domain, &
Ice%fast_domain)
call redistribute_IST_to_IST(Ice%sCS%IST, IST_null, Ice%slow_domain, &
Ice%fast_domain)
else
call SIS_error(FATAL, "Either the pointer to Ice%sCS or the pointer to "//&
"Ice%fCS must be associated in exchange_slow_to_fast_ice.")
endif
else
call SIS_error(FATAL, "exchange_slow_to_fast_ice called with an unrecognized Ice%xtype value.")
endif
call cpu_clock_end(ice_clock_exchange) ; call cpu_clock_end(iceClock)
end subroutine exchange_slow_to_fast_ice
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
!> This subroutine copies information from an ocean_ice_boundary_type into the
!! slow part of an ice_data type, using a coupler-friendly interface.
subroutine unpack_ocean_ice_boundary(Ocean_boundary, Ice)
type(ocean_ice_boundary_type), &
intent(inout) :: Ocean_boundary !< A structure containing information about
!! the ocean that is being shared with the sea-ice.
type(ice_data_type), &
intent(inout) :: Ice !< The publicly visible ice data type in the slow part
!! of which the ocean surface information is to be stored.
if (.not.associated(Ice%sCS)) call SIS_error(FATAL, &
"The pointer to Ice%sCS must be associated in unpack_ocean_ice_boundary.")
call unpack_ocn_ice_bdry(Ocean_boundary, Ice%sCS%OSS, Ice%sCS%IST%ITV, Ice%sCS%G, Ice%sCS%US, &
Ice%sCS%specified_ice, Ice%ocean_fields)
call translate_OSS_to_sOSS(Ice%sCS%OSS, Ice%sCS%IST, Ice%sCS%sOSS, Ice%sCS%G, Ice%sCS%US)
end subroutine unpack_ocean_ice_boundary
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
!> This subroutine converts the information in a publicly visible
!! ocean_ice_boundary_type into an internally visible ocean_sfc_state_type
!! variable.
subroutine unpack_ocn_ice_bdry(OIB, OSS, ITV, G, US, specified_ice, ocean_fields)
type(ocean_ice_boundary_type), intent(in) :: OIB !< A type containing ocean surface fields that
!! are used to drive the sea ice
type(ocean_sfc_state_type), intent(inout) :: OSS !< A structure containing the arrays that describe
!! the ocean's surface state for the ice model.
type(ice_thermo_type), intent(in) :: ITV !< The ice thermodynamics parameter structure.
type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type
type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors
logical, intent(in) :: specified_ice !< If true, use specified ice properties.
type(coupler_3d_bc_type), intent(inout) :: ocean_fields !< A structure of ocean fields, often
!! related to passive tracers.
real, dimension(G%isd:G%ied, G%jsd:G%jed) :: u_nonsym, v_nonsym ! Nonsymmetric velocities [L T-1 ~> m s-1]
real, parameter :: T_0degC = 273.15 ! 0 degrees C in Kelvin
logical :: Cgrid_ocn
integer :: i, j, k, m, n, i2, j2, k2, isc, iec, jsc, jec, i_off, j_off, index
integer :: isd, ied, jsd, jed
isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec
isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
i_off = LBOUND(OIB%t,1) - G%isc ; j_off = LBOUND(OIB%t,2) - G%jsc
call cpu_clock_begin(iceClock) ; call cpu_clock_begin(ice_clock_slow)
! Pass the ocean state through ice on partition 0, unless using specified ice.
if (.not. specified_ice) then
!$OMP parallel do default(shared) private(i2,j2)
do j=jsc,jec ; do i=isc,iec ; i2 = i+i_off ; j2 = j+j_off
OSS%SST_C(i,j) = US%degC_to_C * (OIB%t(i2,j2) - T_0degC)
enddo ; enddo
endif
!$OMP parallel do default(shared) private(i2,j2)
do j=jsc,jec ; do i=isc,iec ; i2 = i+i_off ; j2 = j+j_off
OSS%s_surf(i,j) = US%ppt_to_S*OIB%s(i2,j2)
OSS%T_fr_ocn(i,j) = T_Freeze(OSS%s_surf(i,j), ITV)
OSS%bheat(i,j) = OSS%kmelt*(OSS%SST_C(i,j) - OSS%T_fr_ocn(i,j))
OSS%frazil(i,j) = US%W_m2_to_QRZ_T*US%s_to_T*OIB%frazil(i2,j2)
OSS%sea_lev(i,j) = US%m_to_Z*OIB%sea_level(i2,j2)
enddo ; enddo
Cgrid_ocn = (allocated(OSS%u_ocn_C) .and. allocated(OSS%v_ocn_C))
! Unpack the ocean surface velocities.
if (OIB%stagger == AGRID) then
u_nonsym(:,:) = 0.0 ; v_nonsym(:,:) = 0.0
do j=jsc,jec ; do i=isc,iec ; i2 = i+i_off ; j2 = j+j_off
u_nonsym(i,j) = US%m_s_to_L_T*OIB%u(i2,j2) ; v_nonsym(i,j) = US%m_s_to_L_T*OIB%v(i2,j2)
enddo ; enddo
call pass_vector(u_nonsym, v_nonsym, G%Domain_aux, stagger=AGRID)
if (Cgrid_ocn) then
do j=jsc,jec ; do I=isc-1,iec
OSS%u_ocn_C(I,j) = 0.5*(u_nonsym(i,j) + u_nonsym(i+1,j))
enddo ; enddo
do J=jsc-1,jec ; do i=isc,iec
OSS%v_ocn_C(i,J) = 0.5*(v_nonsym(i,j) + v_nonsym(i,j+1))
enddo ; enddo
else
do J=jsc-1,jec ; do I=isc-1,iec
OSS%u_ocn_B(I,J) = 0.25*((u_nonsym(i,j) + u_nonsym(i+1,j+1)) + &
(u_nonsym(i+1,j) + u_nonsym(i,j+1)))
OSS%v_ocn_B(I,J) = 0.25*((v_nonsym(i,j) + v_nonsym(i+1,j+1)) + &
(v_nonsym(i+1,j) + v_nonsym(i,j+1)))
enddo ; enddo
endif
elseif (OIB%stagger == BGRID_NE) then
if (Cgrid_ocn) then
u_nonsym(:,:) = 0.0 ; v_nonsym(:,:) = 0.0
do j=jsc,jec ; do i=isc,iec ; i2 = i+i_off ; j2 = j+j_off
u_nonsym(i,j) = US%m_s_to_L_T*OIB%u(i2,j2) ; v_nonsym(i,j) = US%m_s_to_L_T*OIB%v(i2,j2)
enddo ; enddo
call pass_vector(u_nonsym, v_nonsym, G%Domain_aux, stagger=BGRID_NE)
do j=jsc,jec ; do I=isc-1,iec
OSS%u_ocn_C(I,j) = 0.5*(u_nonsym(I,J) + u_nonsym(I,J-1))
enddo ; enddo
do J=jsc-1,jec ; do i=isc,iec
OSS%v_ocn_C(i,J) = 0.5*(v_nonsym(I,J) + v_nonsym(I-1,J))
enddo ; enddo
else
do J=jsc,jec ; do I=isc,iec ; i2 = i+i_off ; j2 = j+j_off
OSS%u_ocn_B(I,J) = US%m_s_to_L_T*OIB%u(i2,j2)
OSS%v_ocn_B(I,J) = US%m_s_to_L_T*OIB%v(i2,j2)
enddo ; enddo
if (G%symmetric) &
call fill_symmetric_edges(OSS%u_ocn_B, OSS%v_ocn_B, G%Domain, stagger=BGRID_NE)
endif
elseif (OIB%stagger == CGRID_NE) then
if (Cgrid_ocn) then
do j=jsc,jec ; do I=isc,iec ; i2 = i+i_off ; j2 = j+j_off
OSS%u_ocn_C(I,j) = US%m_s_to_L_T*OIB%u(i2,j2)
enddo ; enddo
do J=jsc,jec ; do i=isc,iec ; i2 = i+i_off ; j2 = j+j_off
OSS%v_ocn_C(i,J) = US%m_s_to_L_T*OIB%v(i2,j2)
enddo ; enddo
if (G%symmetric) &
call fill_symmetric_edges(OSS%u_ocn_C, OSS%v_ocn_C, G%Domain, stagger=CGRID_NE)
else
u_nonsym(:,:) = 0.0 ; v_nonsym(:,:) = 0.0
do j=jsc,jec ; do i=isc,iec ; i2 = i+i_off ; j2 = j+j_off
u_nonsym(I,j) = US%m_s_to_L_T*OIB%u(i2,j2) ; v_nonsym(i,J) = US%m_s_to_L_T*OIB%v(i2,j2)
enddo ; enddo
call pass_vector(u_nonsym, v_nonsym, G%Domain_aux, stagger=CGRID_NE)
do J=jsc-1,jec ; do I=isc-1,iec
OSS%u_ocn_B(I,J) = 0.5*(u_nonsym(I,j) + u_nonsym(I,j+1))
OSS%v_ocn_B(I,J) = 0.5*(v_nonsym(i,J) + v_nonsym(i+1,J))
enddo ; enddo
endif
else
call SIS_error(FATAL, "unpack_ocn_ice_bdry: Unrecognized OIB%stagger.")
endif
! Fill in the halo values.
if (Cgrid_ocn) then
call pass_vector(OSS%u_ocn_C, OSS%v_ocn_C, G%Domain, stagger=CGRID_NE)
else
call pass_vector(OSS%u_ocn_B, OSS%v_ocn_B, G%Domain, stagger=BGRID_NE)
endif
call pass_var(OSS%sea_lev, G%Domain)
! Transfer the ocean state for extra tracer fluxes.
call coupler_type_spawn(OIB%fields, OSS%tr_fields, (/isd, isc, iec, ied/), &
(/jsd, jsc, jec, jed/), as_needed=.true. )
call coupler_type_copy_data(OIB%fields, OSS%tr_fields)
call cpu_clock_end(ice_clock_slow) ; call cpu_clock_end(iceClock)
end subroutine unpack_ocn_ice_bdry