forked from NOAA-GFDL/SIS2
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathice_model.F90
2587 lines (2178 loc) · 119 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
!***********************************************************************
!* 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 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 ([email protected]). !
! SIS2 has been developed by Robert Hallberg and Mike Winton, with !
! contributions from many people at NOAA/GFDL, including Alistair Adcroft and !
! Niki Zadeh. !
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
module ice_model_mod
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_diag_mediator, only : query_SIS_averaging_enabled, SIS_diag_ctrl
! use SIS_diag_mediator, only : register_diag_field=>register_SIS_diag_field
use SIS_get_input, only : Get_SIS_input, directories
use SIS_sum_output, only : SIS_sum_output_init, write_ice_statistics
use SIS_transcribe_grid, only : copy_dyngrid_to_SIS_horgrid, copy_SIS_horgrid_to_dyngrid
use MOM_checksums, only : chksum, uchksum, vchksum, Bchksum
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_obsolete_params, only : obsolete_logical
use MOM_string_functions, only : uppercase
use MOM_time_manager, only : time_type, time_type_to_real, real_to_time_type
use MOM_time_manager, only : set_date, set_time, operator(+), operator(-)
use MOM_time_manager, only : operator(>), operator(*), operator(/), operator(/=)
use fms_mod, only : file_exist, clock_flag_default
use fms_io_mod, only : set_domain, nullify_domain, restore_state, query_initialized
use fms_io_mod, only : restore_state, query_initialized
use fms_io_mod, only : register_restart_field, restart_file_type
use mpp_mod, only : mpp_clock_id, mpp_clock_begin, mpp_clock_end
use mpp_mod, only : CLOCK_COMPONENT, CLOCK_ROUTINE
use mpp_domains_mod, only : mpp_broadcast_domain
use astronomy_mod, only : astronomy_init, astronomy_end
use astronomy_mod, only : universal_time, orbital_time, diurnal_solar, daily_mean_solar
use coupler_types_mod, only : coupler_3d_bc_type
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_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 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 SIS_ctrl_types, only : SIS_slow_CS, SIS_fast_CS
use SIS_ctrl_types, only : ice_diagnostics_init, ice_diags_fast_init
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
use SIS_types, only : fast_ice_avg_type, alloc_fast_ice_avg, dealloc_fast_ice_avg
use SIS_types, only : ice_rad_type, ice_rad_register_restarts, dealloc_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 : copy_IST_to_IST, copy_FIA_to_FIA, copy_sOSS_to_sOSS
use SIS_types, only : redistribute_IST_to_IST, redistribute_FIA_to_FIA
use SIS_types, only : redistribute_sOSS_to_sOSS, FIA_chksum, IOF_chksum
use ice_utils_mod, only : post_avg, ice_grid_chksum
use SIS_hor_grid, only : SIS_hor_grid_type, set_hor_grid, SIS_hor_grid_end, set_first_direction
use SIS_fixed_initialization, only : SIS_initialize_fixed
use ice_grid, only : set_ice_grid, ice_grid_end, ice_grid_type
use ice_spec_mod, only : get_sea_surface
use SIS_tracer_registry, only : register_SIS_tracer, register_SIS_tracer_pair
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 ice_thm_mod, only : slab_ice_optics
use SIS_dyn_trans, only : SIS_dynamics_trans, update_icebergs
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_transport_CS, SIS_dyn_trans_sum_output_CS
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_fast_thermo, only : do_update_ice_model_fast, avg_top_quantities
use SIS_fast_thermo, only : SIS_fast_thermo_init, SIS_fast_thermo_end
use SIS_optics, only : ice_optics_SIS2, SIS_optics_init, SIS_optics_end
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, get_SIS2_thermo_coefs
use SIS2_ice_thm, only : enth_from_TS, Temp_from_En_S, T_freeze, ice_thermo_type
use ice_bergs, only : icebergs, icebergs_run, icebergs_init, icebergs_end
implicit none ; private
public :: ice_data_type, ocean_ice_boundary_type, atmos_ice_boundary_type, land_ice_boundary_type
public :: ice_model_init, ice_model_end, update_ice_model_fast, ice_stock_pe
public :: update_ice_model_slow_up, update_ice_model_slow_dn ! The old Verona interfaces.
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 :: unpack_ocean_ice_boundary, exchange_slow_to_fast_ice, set_ice_surface_fields
public :: ice_model_fast_cleanup, unpack_land_ice_boundary
public :: exchange_fast_to_slow_ice, update_ice_model_slow
integer :: iceClock, iceClock1, iceCLock2, iceCLock3
integer, parameter :: REDIST=2, DIRECT=3
contains
!-----------------------------------------------------------------------
!> Update the sea-ice state due to slow processes, including dynamics,
!! freezing and melting, precipitation, and transport.
subroutine update_ice_model_slow_dn ( Atmos_boundary, Land_boundary, Ice )
type(atmos_ice_boundary_type), &
intent(in) :: Atmos_boundary !< Atmos_boundary is not actually used, and
!! is still here only for backward compatibilty with the
!! interface to Verona and earlier couplers.
type(land_ice_boundary_type), &
intent(in) :: Land_boundary !< A structure containing information about
!! the fluxes from the land that is being shared with the
!! sea-ice. If this argument is not present, it is assumed
!! that this information has already been exchanged.
type(ice_data_type), &
intent(inout) :: Ice !< The publicly visible ice data type; this must always be
!! present, but is optional because of an unfortunate
!! order of arguments.
if (.not.associated(Ice%sCS)) call SIS_error(FATAL, &
"The pointer to Ice%sCS must be associated in update_ice_model_slow_dn.")
call mpp_clock_begin(iceClock) ; call mpp_clock_begin(iceClock2)
call ice_model_fast_cleanup(Ice)
call unpack_land_ice_boundary(Ice, Land_boundary)
! In the case where fast and slow ice PEs are not the same, this call would
! need to be replaced by a routine that does inter-processor receives.
call exchange_fast_to_slow_ice(Ice)
call mpp_clock_end(iceClock2) ; call mpp_clock_end(iceClock)
call update_ice_model_slow(Ice)
end subroutine update_ice_model_slow_dn
!-----------------------------------------------------------------------
!> 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.
real :: dt_slow ! The time step over which to advance the model.
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_model_slow.")
call mpp_clock_begin(iceClock) ; call mpp_clock_begin(iceClock2)
! 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 = time_type_to_real(Ice%sCS%Time_step_slow)
if (Ice%sCS%debug) then
call Ice_public_type_chksum("Start update_ice_model_slow", Ice)
call FIA_chksum("Start update_ice_model_slow", Ice%sCS%FIA, Ice%sCS%G)
call IOF_chksum("Start update_ice_model_slow", Ice%sCS%IOF, Ice%sCS%G)
endif
! Store some diagnostic fluxes...
!$OMP parallel do default(none) shared(Ice)
do j=Ice%sCS%G%jsc,Ice%sCS%G%jec ; do i=Ice%sCS%G%isc,Ice%sCS%G%iec
Ice%sCS%FIA%calving_preberg(i,j) = Ice%sCS%FIA%calving(i,j)
Ice%sCS%FIA%calving_hflx_preberg(i,j) = Ice%sCS%FIA%calving_hflx(i,j)
enddo ; enddo
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) - Ice%sCS%G%isc
j_off = LBOUND(Ice%flux_t,2) - Ice%sCS%G%jsc
!$OMP parallel do default(none) shared(Ice,i_off,j_off) private(i2,j2)
do j=Ice%sCS%G%jsc,Ice%sCS%G%jec ; do i=Ice%sCS%G%isc,Ice%sCS%G%iec
i2 = i+i_off ; j2 = j+j_off
Ice%sCS%IOF%flux_u_ocn(i,j) = Ice%flux_u(i2,j2)
Ice%sCS%IOF%flux_v_ocn(i,j) = Ice%flux_v(i2,j2)
enddo ; enddo
endif
call mpp_clock_end(iceClock2) ; call mpp_clock_end(iceClock)
call update_icebergs(Ice%sCS%IST, Ice%sCS%OSS, Ice%sCS%IOF, Ice%sCS%FIA, Ice%icebergs, &
dt_slow, Ice%sCS%G, Ice%sCS%IG, Ice%sCS%dyn_trans_CSp)
call mpp_clock_begin(iceClock) ; call mpp_clock_begin(iceClock2)
if (Ice%sCS%debug) then
call FIA_chksum("After update_icebergs", Ice%sCS%FIA, Ice%sCS%G)
endif
endif
if (Ice%sCS%debug) then
call Ice_public_type_chksum("Before slow_thermodynamics", Ice)
call IOF_chksum("Before slow_thermodynamics", Ice%sCS%IOF, Ice%sCS%G)
endif
call slow_thermodynamics(Ice%sCS%IST, dt_slow, Ice%sCS%slow_thermo_CSp, &
Ice%sCS%OSS, Ice%sCS%FIA, Ice%sCS%IOF, Ice%sCS%G, Ice%sCS%IG)
! 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.
call pass_vector(Ice%sCS%FIA%WindStr_x, Ice%sCS%FIA%WindStr_y, &
Ice%sCS%G%Domain, stagger=AGRID, complete=.false.)
call pass_vector(Ice%sCS%FIA%WindStr_ocn_x, Ice%sCS%FIA%WindStr_ocn_y, &
Ice%sCS%G%Domain, stagger=AGRID)
call pass_var(Ice%sCS%FIA%ice_cover, Ice%sCS%G%Domain, complete=.false.)
call pass_var(Ice%sCS%FIA%ice_free, Ice%sCS%G%Domain, complete=.true.)
call pass_var(Ice%sCS%IST%part_size, Ice%sCS%G%Domain)
call pass_var(Ice%sCS%IST%mH_ice, Ice%sCS%G%Domain, complete=.false.)
call pass_var(Ice%sCS%IST%mH_pond, Ice%sCS%G%Domain, complete=.false.)
call pass_var(Ice%sCS%IST%mH_snow, Ice%sCS%G%Domain, complete=.true.)
if (Ice%sCS%debug) then
call Ice_public_type_chksum("Before SIS_dynamics_trans", Ice)
call IOF_chksum("Before SIS_dynamics_trans", Ice%sCS%IOF, Ice%sCS%G)
endif
call SIS_dynamics_trans(Ice%sCS%IST, Ice%sCS%OSS, Ice%sCS%FIA, Ice%sCS%IOF, &
dt_slow, Ice%sCS%dyn_trans_CSp, Ice%icebergs, Ice%sCS%G, Ice%sCS%IG)
if (Ice%sCS%debug) then
call Ice_public_type_chksum("Before set_ocean_top_fluxes", Ice)
call IOF_chksum("Before set_ocean_top_fluxes", Ice%sCS%IOF, Ice%sCS%G)
call IST_chksum("Before set_ocean_top_fluxes", Ice%sCS%IST, Ice%sCS%G, Ice%sCS%IG)
endif
! Set up the thermodynamic fluxes in the externally visible structure Ice.
call set_ocean_top_fluxes(Ice, Ice%sCS%IST, Ice%sCS%IOF, Ice%sCS%FIA, Ice%sCS%OSS, &
Ice%sCS%G, Ice%sCS%IG, Ice%sCS)
if (Ice%sCS%debug) then
call Ice_public_type_chksum("End update_ice_model_slow", Ice)
endif
!### THIS NO LONGER WORKS ON SLOW ICE PES.
! if (Ice%sCS%bounds_check) then
! call Ice_public_type_bounds_check(Ice, Ice%sCS%G, "End update_ice_slow")
! endif
call mpp_clock_end(iceClock2) ; call mpp_clock_end(iceClock)
end subroutine update_ice_model_slow
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
!> 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%part_size, &
Ice%fCS%G, Ice%fCS%IG)
call store_Tskin(Ice%fCS%IST, Ice%fCS%sOSS, Ice%fCS%Rad, Ice%fCS%G, Ice%fCS%IG)
end subroutine ice_model_fast_cleanup
!> store_Tskin stores a filled-in array of actual, interpolated or plausible
!> ice/snow-surface skin temperatures for later use.
subroutine store_Tskin(IST, OSS, Rad, G, IG)
type(ice_state_type), intent(in) :: IST !< The ice state type whose values of t_surf are being used.
type(simple_OSS_type), intent(in) :: OSS !< An ocean surface state type, used for the surface freezing point.
type(ice_rad_type), intent(inout) :: Rad !< An ice rad type, where the t_skin are being stored.
type(SIS_hor_grid_type), intent(in) :: G !< The sea-ice lateral grid type.
type(ice_grid_type), intent(in) :: IG !< The sea-ice grid type.
integer, dimension(G%isd:G%ied) :: k_thick, k_thin
real, dimension(G%isd:G%ied) :: &
t_thick, t_thin, mH_thin, t_fr_ocn
real :: wt2, mH_cat_tgt
real, parameter :: T_0degC = 273.15 ! 0 degrees C in Kelvin
logical, dimension(G%isd:G%ied) :: infill
logical :: any_ice
integer :: i, j, k, k2, k3, isc, iec, jsc, jec, ncat, NkIce
isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec
ncat = IG%CatIce ; NkIce = IG%NkIce
do j=jsc,jec
! Determine which categories exist.
do i=isc,iec ; k_thick(i) = -1 ; k_thin(i) = ncat+1 ; mH_thin(i) = 0.0 ; enddo
any_ice = .false.
do k=1,ncat ; do i=isc,iec ; if (IST%part_size(i,j,k) > 0.0) then
k_thick(i) = k
t_thick(i) = IST%t_surf(i,j,k)
any_ice = .true.
endif ; enddo ; enddo
do i=isc,iec
T_fr_ocn(i) = G%mask2dT(i,j) * T_Freeze(OSS%s_surf(i,j), IST%ITV)
if (k_thick(i) < 1) t_thick(i) = T_fr_ocn(i) + T_0degC
enddo
if (.not.any_ice) then
! This entire row is ice free.
do k=1,ncat ; do i=isc,iec
Rad%t_skin(i,j,k) = G%mask2dT(i,j)*T_fr_ocn(i) + T_0degC
enddo ; enddo
else
do k=ncat,1,-1 ; do i=isc,iec ; if (IST%part_size(i,j,k) > 0.0) then
k_thin(i) = k
t_thin(i) = IST%t_surf(i,j,k)
mH_thin(i) = IST%mH_ice(i,j,k)
endif ; enddo ; enddo
! The logic for the k=1 (thinest) category is particularly simple.
do i=isc,iec ; if (G%mask2dT(i,j) > 0.0) then
if (IST%part_size(i,j,1) > 0.0) then
Rad%t_skin(i,j,1) = IST%t_surf(i,j,1)
else
Rad%t_skin(i,j,1) = T_fr_ocn(i) + T_0degC
endif
endif ; enddo
do k=2,ncat ; do i=isc,iec
if (G%mask2dT(i,j) > 0.0) then
if (k > k_thick(i)) then
! This is the most common case, since it applies to ice-free water.
Rad%t_skin(i,j,k) = t_thick(i)
elseif (IST%part_size(i,j,k) > 0.0) then
! Copy over the valid value - no interpolation is needed.
Rad%t_skin(i,j,k) = IST%t_surf(i,j,k)
elseif (k < k_thin(i)) then
! Linearly interpolate to the ocean's freezing temperature.
mH_cat_tgt = 0.5*(IG%mH_cat_bound(K) + IG%mH_cat_bound(K+1))
wt2 = 0.0 ! If in doubt, use t_thin. This should not happen.
if (mH_thin(i) > mH_cat_tgt) wt2 = (mH_thin(i) - mH_cat_tgt) / mH_thin(i)
Rad%t_skin(i,j,k) = wt2*(T_fr_ocn(i) + T_0degC) + (1.0-wt2)*t_thin(i)
else
! Linearly interpolate between bracketing neighboring
! categories with valid values.
do k2=k+1,ncat ; if (IST%part_size(i,j,k2) > 0.0) exit ; enddo
do k3=k-1,1,-1 ; if (IST%part_size(i,j,k3) > 0.0) exit ; enddo
mH_cat_tgt = 0.5*(IG%mH_cat_bound(K) + IG%mH_cat_bound(K+1))
wt2 = 0.5
if ((IST%mH_ice(i,j,k3) > mH_cat_tgt) .and. &
(IST%mH_ice(i,j,k2) < mH_cat_tgt)) &
wt2 = (IST%mH_ice(i,j,k3) - mH_cat_tgt) / &
(IST%mH_ice(i,j,k3) - IST%mH_ice(i,j,k2))
Rad%t_skin(i,j,k) = wt2*IST%t_surf(i,j,k2) + (1.0-wt2) * IST%t_surf(i,j,k3)
endif
else ! This is a land point.
Rad%t_skin(i,j,k) = T_0degC
endif
enddo ; enddo
endif
enddo
end subroutine store_Tskin
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
!> 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()
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
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) &
!$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) = LIB%runoff(i2,j2)
FIA%calving(i,j) = LIB%calving(i2,j2)
FIA%runoff_hflx(i,j) = LIB%runoff_hflx(i2,j2)
FIA%calving_hflx(i,j) = 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)
endif
end subroutine unpack_land_ice_boundary
!> This subroutine copies information (mostly fluxes and the updated tempertures)
!! 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.
if (.not.associated(Ice%fCS) .and. .not.associated(Ice%sCS)) call SIS_error(FATAL, &
"For now, 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 SIS_mesg("Copying Ice%fCS%FIA to Ice%sCS%FIA in update_ice_model_slow_dn.")
if(Ice%xtype == DIRECT) then
call copy_FIA_to_FIA(Ice%fCS%FIA, Ice%sCS%FIA, Ice%fCS%G%HI, Ice%sCS%G%HI, &
Ice%sCS%IG)
else if(Ice%xtype == REDIST) then
call redistribute_FIA_to_FIA(Ice%fCS%FIA, Ice%sCS%FIA, Ice%fast_domain, Ice%slow_domain, &
Ice%sCS%G, Ice%sCS%IG)
endif
endif
if (.not.associated(Ice%fCS%IST, Ice%sCS%IST)) then
! call SIS_mesg("Copying Ice%fCS%IST to Ice%sCS%IST in update_ice_model_slow_dn.")
if(Ice%xtype == DIRECT) then
call copy_IST_to_IST(Ice%fCS%IST, Ice%sCS%IST, Ice%fCS%G%HI, Ice%sCS%G%HI, &
Ice%fCS%IG)
else if(Ice%xtype == REDIST) then
call redistribute_IST_to_IST(Ice%fCS%IST, Ice%sCS%IST, Ice%fast_domain, &
Ice%slow_domain)
endif
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, IG, sCS)
type(ice_data_type), intent(inout) :: Ice
type(ice_state_type), intent(inout) :: IST
type(ice_ocean_flux_type), intent(in) :: IOF
type(fast_ice_avg_type), intent(in) :: FIA
type(ocean_sfc_state_type), intent(in) :: OSS
type(SIS_hor_grid_type), intent(inout) :: G
type(ice_grid_type), intent(in) :: IG
type(SIS_slow_CS), intent(in) :: sCS
real :: I_count
integer :: i, j, k, isc, iec, jsc, jec, m, n
integer :: i2, j2, i_off, j_off, ind, 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)
call IOF_chksum("Start set_ocean_top_fluxes", IOF, G)
call FIA_chksum("Start set_ocean_top_fluxes", FIA, G)
endif
! This block of code is probably unneccessary.
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
do n=1,Ice%ocean_fluxes%num_bcs ; do m=1,Ice%ocean_fluxes%bc(n)%num_fields
Ice%ocean_fluxes%bc(n)%field(m)%values(:,:) = 0.0
enddo ; enddo
! 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(none) shared(isc,iec,jsc,jec,ncat,Ice,IST,IG,i_off,j_off) &
!$OMP 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) * &
(IG%H_to_kg_m2 * (IST%mH_snow(i,j,k) + IST%mH_ice(i,j,k)))
enddo ; 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(none) shared(isc,iec,jsc,jec,Ice,IST,IOF,FIA,i_off,j_off,G,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_u(i2,j2) = IOF%flux_u_ocn(i,j)
Ice%flux_v(i2,j2) = IOF%flux_v_ocn(i,j)
Ice%flux_t(i2,j2) = IOF%flux_t_ocn_top(i,j)
Ice%flux_q(i2,j2) = IOF%flux_q_ocn_top(i,j)
Ice%flux_sw_vis_dir(i2,j2) = IOF%flux_sw_vis_dir_ocn(i,j)
Ice%flux_sw_vis_dif(i2,j2) = IOF%flux_sw_vis_dif_ocn(i,j)
Ice%flux_sw_nir_dir(i2,j2) = IOF%flux_sw_nir_dir_ocn(i,j)
Ice%flux_sw_nir_dif(i2,j2) = IOF%flux_sw_nir_dif_ocn(i,j)
Ice%flux_lw(i2,j2) = IOF%flux_lw_ocn_top(i,j)
Ice%flux_lh(i2,j2) = IOF%flux_lh_ocn_top(i,j)
Ice%fprec(i2,j2) = IOF%fprec_ocn_top(i,j)
Ice%lprec(i2,j2) = IOF%lprec_ocn_top(i,j)
Ice%runoff(i2,j2) = FIA%runoff(i,j)
Ice%calving(i2,j2) = FIA%calving(i,j)
Ice%runoff_hflx(i2,j2) = FIA%runoff_hflx(i,j)
Ice%calving_hflx(i2,j2) = FIA%calving_hflx(i,j)
Ice%flux_salt(i2,j2) = IOF%flux_salt(i,j)
Ice%SST_C(i2,j2) = OSS%t_ocn(i,j)
if (IOF%slp2ocean) then
Ice%p_surf(i2,j2) = FIA%p_atm_surf(i,j) - 1e5 ! SLP - 1 std. atmosphere, in Pa.
else
Ice%p_surf(i2,j2) = 0.0
endif
Ice%p_surf(i2,j2) = Ice%p_surf(i2,j2) + 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) + IOF%melt_nudge(i,j)
enddo ; enddo
endif
ind = 0
do n=1,Ice%ocean_fluxes%num_bcs ; do m=1,Ice%ocean_fluxes%bc(n)%num_fields
ind = ind + 1
do j=jsc,jec ; do i=isc,iec
i2 = i+i_off ; j2 = j+j_off ! Use these to correct for indexing differences.
Ice%ocean_fluxes%bc(n)%field(m)%values(i2,j2) = IOF%tr_flux_ocn_top(i,j,ind)
enddo ; enddo
enddo ; enddo
! This extra block is required with the Verona and earlier versions of the coupler.
i_off = LBOUND(Ice%part_size,1) - G%isc ; j_off = LBOUND(Ice%part_size,2) - G%jsc
if (Ice%shared_slow_fast_PEs) then
if ((Ice%fCS%G%iec-Ice%fCS%G%isc==iec-isc) .and. &
(Ice%fCS%G%jec-Ice%fCS%G%jsc==jec-jsc)) then
! The fast and slow ice PEs are using the same PEs and layout, so the
! part_size arrays can be copied directly from the fast ice PEs.
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,ncat,Ice,IST,i_off,j_off) &
!$OMP private(i2,j2)
do j=jsc,jec ; do k=0,ncat ; do i=isc,iec
i2 = i+i_off ; j2 = j+j_off! Use these to correct for indexing differences.
Ice%part_size(i2,j2,k+1) = IST%part_size(i,j,k)
enddo ; enddo ; enddo
endif
endif
if (sCS%debug) then
call Ice_public_type_chksum("End set_ocean_top_fluxes", Ice)
call IST_chksum("Start set_ocean_top_fluxes", IST, G, IG)
endif
end subroutine set_ocean_top_fluxes
! Coupler interface to provide ocean surface data to atmosphere.
!
!> update_ice_model_slow_up prepares the ice surface data for forcing the atmosphere
!! and also unpacks the data from the ocean and shares it between the fast and
!! slow ice structures.
subroutine update_ice_model_slow_up ( 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.
if (.not.associated(Ice%fCS)) call SIS_error(FATAL, &
"The pointer to Ice%fCS must be associated in update_ice_model_slow_up.")
if (.not.associated(Ice%sCS)) call SIS_error(FATAL, &
"The pointer to Ice%sCS must be associated in update_ice_model_slow_up.")
call unpack_ocn_ice_bdry(Ocean_boundary, Ice%sCS%OSS, Ice%sCS%G, &
Ice%sCS%IST%t_surf(:,:,0), 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%IST%ITV)
call exchange_slow_to_fast_ice(Ice)
call set_ice_surface_fields(Ice)
end subroutine update_ice_model_slow_up
!> 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.
call mpp_clock_begin(iceClock) ; call mpp_clock_begin(iceClock1)
if (.not.associated(Ice%fCS) .and. .not.associated(Ice%sCS)) call SIS_error(FATAL, &
"For now, 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
if(Ice%xtype == DIRECT) then
call copy_sOSS_to_sOSS(Ice%sCS%sOSS, Ice%fCS%sOSS, Ice%sCS%G%HI, Ice%fCS%G%HI)
else if(Ice%xtype == REDIST) then
call redistribute_sOSS_to_sOSS(Ice%sCS%sOSS, Ice%fCS%sOSS, Ice%slow_domain, &
Ice%fast_domain, Ice%fCS%G%HI)
endif
endif
if (.not.associated(Ice%fCS%IST, Ice%sCS%IST)) then
! call SIS_mesg("Copying Ice%sCS%IST to Ice%fCS%IST in update_ice_model_slow_up.")
if(Ice%xtype == DIRECT) then
call copy_IST_to_IST(Ice%sCS%IST, Ice%fCS%IST, Ice%sCS%G%HI, Ice%fCS%G%HI, &
Ice%sCS%IG)
else
call redistribute_IST_to_IST(Ice%sCS%IST, Ice%fCS%IST, Ice%slow_domain, &
Ice%fast_domain)
endif
endif
call mpp_clock_end(iceClock1) ; call mpp_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%G, &
Ice%sCS%IST%t_surf(:,:,0), 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%IST%ITV)
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, G, t_surf_ocn_K, specified_ice, ocean_fields)
type(ocean_ice_boundary_type), intent(in) :: OIB
type(ocean_sfc_state_type), intent(inout) :: OSS
type(SIS_hor_grid_type), intent(inout) :: G
real, dimension(G%isd:G%ied, G%jsd:G%jed), &
intent(inout) :: t_surf_ocn_K ! The ocean surface temperature in Kelvin.
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
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
isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec
i_off = LBOUND(OIB%t,1) - G%isc ; j_off = LBOUND(OIB%t,2) - G%jsc
call mpp_clock_begin(iceClock) ; call mpp_clock_begin(iceClock1)
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,OSS,OIB,i_off,j_off) &
!$OMP private(i2,j2)
do j=jsc,jec ; do i=isc,iec ; i2 = i+i_off ; j2 = j+j_off
OSS%t_ocn(i,j) = OIB%t(i2,j2) - T_0degC
OSS%s_surf(i,j) = OIB%s(i2,j2)
OSS%frazil(i,j) = OIB%frazil(i2,j2)
OSS%sea_lev(i,j) = OIB%sea_level(i2,j2)
enddo ; enddo
! Pass the ocean state through ice on partition 0, unless using specified ice.
if (.not. specified_ice) then
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,t_surf_ocn_K,OIB,i_off,j_off) &
!$OMP private(i2,j2)
do j=jsc,jec ; do i=isc,iec ; i2 = i+i_off ; j2 = j+j_off
t_surf_ocn_K(i,j) = OIB%t(i2,j2)
enddo ; enddo
endif
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) = OIB%u(i2,j2) ; v_nonsym(i,j) = 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) = OIB%u(i2,j2) ; v_nonsym(i,j) = 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) = OIB%u(i2,j2)
OSS%v_ocn_B(I,J) = 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) = 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) = 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) = OIB%u(i2,j2) ; v_nonsym(i,J) = 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.
if (OSS%num_tr<0) then
! Determine the number of tracer arrays and allocate them.
OSS%num_tr = 0
do n=1,OIB%fields%num_bcs
OSS%num_tr = OSS%num_tr + OIB%fields%bc(n)%num_fields
enddo
if (OSS%num_tr > 0) then
allocate(OSS%tr_array(G%isd:G%ied,G%jsd:G%jed,OSS%num_tr)) ; OSS%tr_array(:,:,:) = 0.0
endif
endif
index = 0
do n=1,OIB%fields%num_bcs ; do m=1,OIB%fields%bc(n)%num_fields
index = index+1
if (index == 1) then
i_off = LBOUND(OIB%fields%bc(n)%field(m)%values,1) - G%isc
j_off = LBOUND(OIB%fields%bc(n)%field(m)%values,2) - G%jsc
endif
do j=jsc,jec ; do i=isc,iec ; i2 = i+i_off ; j2 = j+j_off
OSS%tr_array(i,j,index) = OIB%fields%bc(n)%field(m)%values(i2,j2)
enddo ; enddo
enddo ; enddo
call mpp_clock_end(iceClock1) ; call mpp_clock_end(iceClock)
end subroutine unpack_ocn_ice_bdry
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
!> translate_OSS_to_sOSS translates the full ocean surface state, as seen by the slow
!! ice processors into a simplified version with the fields that are shared with
!! the atmosphere and the fast ice thermodynamics.
subroutine translate_OSS_to_sOSS(OSS, IST, sOSS, G, ITV)
type(ocean_sfc_state_type), intent(in) :: OSS
type(ice_state_type), intent(in) :: IST
type(simple_OSS_type), intent(inout) :: sOSS
type(SIS_hor_grid_type), intent(in) :: G
type(ice_thermo_type), intent(in) :: ITV
integer :: i, j, k, m, n, i2, j2, k2, isc, iec, jsc, jec, i_off, j_off
isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,G,sOSS,OSS,IST,ITV)
do j=jsc,jec ; do i=isc,iec
sOSS%s_surf(i,j) = OSS%s_surf(i,j)
sOSS%t_ocn(i,j) = OSS%t_ocn(i,j)
if (G%mask2dT(i,j) > 0.5) then
sOSS%bheat(i,j) = OSS%kmelt*(OSS%t_ocn(i,j) - T_Freeze(OSS%s_surf(i,j), ITV))
! Interpolate the ocean and ice velocities onto tracer cells.
if (OSS%Cgrid_dyn) then
sOSS%u_ocn_A(i,j) = 0.5*(OSS%u_ocn_C(I,j) + OSS%u_ocn_C(I-1,j))
sOSS%v_ocn_A(i,j) = 0.5*(OSS%v_ocn_C(i,J) + OSS%v_ocn_C(i,J-1))
else
sOSS%u_ocn_A(i,j) = 0.25*((OSS%u_ocn_B(I,J) + OSS%u_ocn_B(I-1,J-1)) + &
(OSS%u_ocn_B(I,J-1) + OSS%u_ocn_B(I-1,J)) )
sOSS%v_ocn_A(i,j) = 0.25*((OSS%v_ocn_B(I,J) + OSS%v_ocn_B(I-1,J-1)) + &
(OSS%v_ocn_B(I,J-1) + OSS%v_ocn_B(I-1,J)) )
endif
if (IST%Cgrid_dyn) then
sOSS%u_ice_A(i,j) = 0.5*(IST%u_ice_C(I,j) + IST%u_ice_C(I-1,j))
sOSS%v_ice_A(i,j) = 0.5*(IST%v_ice_C(i,J) + IST%v_ice_C(i,J-1))
else
sOSS%u_ice_A(i,j) = 0.25*((IST%u_ice_B(I,J) + IST%u_ice_B(I-1,J-1)) + &
(IST%u_ice_B(I,J-1) + IST%u_ice_B(I-1,J)) )
sOSS%v_ice_A(i,j) = 0.25*((IST%v_ice_B(I,J) + IST%v_ice_B(I-1,J-1)) + &
(IST%v_ice_B(I,J-1) + IST%v_ice_B(I-1,J)) )
endif
else ! This is a land point.
sOSS%bheat(i,j) = 0.0
sOSS%u_ocn_A(i,j) = 0.0 ; sOSS%v_ocn_A(i,j) = 0.0
sOSS%u_ice_A(i,j) = 0.0 ; sOSS%v_ice_A(i,j) = 0.0
endif
enddo ; enddo
if (sOSS%num_tr<0) then
sOSS%num_tr = OSS%num_tr
if (sOSS%num_tr > 0) then
allocate(sOSS%tr_array(G%isd:G%ied,G%jsd:G%jed,sOSS%num_tr)) ; sOSS%tr_array(:,:,:) = 0.0
endif
endif
do m=1,OSS%num_tr ; do j=jsc,jec ; do i=isc,iec
sOSS%tr_array(i,j,m) = OSS%tr_array(i,j,m)
enddo ; enddo ; enddo
end subroutine translate_OSS_to_sOSS
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
!> set_ice_surface_fields prepares the ice surface state for atmosphere fast
!! physics and does precalculation of ice radiative properties.
subroutine set_ice_surface_fields(Ice)
type(ice_data_type), intent(inout) :: Ice !< The publicly visible ice data type whose
!! surface properties are being set.
call mpp_clock_begin(iceClock) ; call mpp_clock_begin(iceClock1)
if (.not.associated(Ice%fCS)) call SIS_error(FATAL, &
"The pointer to Ice%fCS must be associated in set_ice_surface_fields.")
call set_ice_surface_state(Ice, Ice%fCS%IST, Ice%fCS%sOSS, Ice%fCS%Rad, &
Ice%fCS%FIA, Ice%fCS%G, Ice%fCS%IG, Ice%fCS )
call mpp_clock_end(iceClock1) ; call mpp_clock_end(iceClock)
end subroutine set_ice_surface_fields
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
! set_ice_surface_state - prepare surface state for atmosphere fast physics !
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
subroutine set_ice_surface_state(Ice, IST, OSS, Rad, FIA, G, IG, fCS)
type(ice_data_type), intent(inout) :: Ice
type(ice_state_type), intent(inout) :: IST
type(simple_OSS_type), intent(in) :: OSS
type(ice_rad_type), intent(inout) :: Rad
type(fast_ice_avg_type), intent(inout) :: FIA
type(SIS_hor_grid_type), intent(inout) :: G
type(ice_grid_type), intent(in) :: IG
type(SIS_fast_CS), intent(inout) :: fCS
real, dimension(G%isc:G%iec,G%jsc:G%jec) :: m_ice_tot
real, dimension(IG%NkIce) :: sw_abs_lay
real :: u, v
real :: area_pt
real :: rho_ice ! The nominal density of sea ice in kg m-3.
real :: rho_snow ! The nominal density of snow in kg m-3.
type(time_type) :: dt_r ! A temporary radiation timestep.
integer :: i, j, k, m, n, i2, j2, k2, isc, iec, jsc, jec, ncat, i_off, j_off
integer :: index
real :: H_to_m_ice ! The specific volumes of ice and snow times the
real :: H_to_m_snow ! conversion factor from thickness units, in m H-1.
real, parameter :: T_0degC = 273.15 ! 0 degrees C in Kelvin
logical :: slab_ice ! If true, use the very old slab ice thermodynamics,
! with effectively zero heat capacity of ice and snow.
isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; ncat = IG%CatIce
i_off = LBOUND(Ice%t_surf,1) - G%isc ; j_off = LBOUND(Ice%t_surf,2) - G%jsc
call get_SIS2_thermo_coefs(IST%ITV, rho_ice=rho_ice, rho_snow=rho_snow, slab_ice=slab_ice)
H_to_m_snow = IG%H_to_kg_m2 / Rho_snow ; H_to_m_ice = IG%H_to_kg_m2 / Rho_ice
if (fCS%bounds_check) &
call IST_bounds_check(IST, G, IG, "Start of set_ice_surface_state") !, OSS=OSS)
if (fCS%debug) then
call IST_chksum("Start set_ice_surface_state", IST, G, IG)
call Ice_public_type_chksum("Start set_ice_surface_state", Ice)
endif
m_ice_tot(:,:) = 0.0
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,G,IST,OSS,FIA,ncat,m_ice_tot)
do j=jsc,jec
do k=1,ncat ; do i=isc,iec
FIA%tmelt(i,j,k) = 0.0 ; FIA%bmelt(i,j,k) = 0.0
m_ice_tot(i,j) = m_ice_tot(i,j) + IST%mH_ice(i,j,k) * IST%part_size(i,j,k)
enddo ; enddo
do i=isc,iec
if (m_ice_tot(i,j) > 0.0) then
FIA%bheat(i,j) = OSS%bheat(i,J)
else
FIA%bheat(i,j) = 0.0
endif
enddo
enddo
if (fCS%Eulerian_tsurf) then
do k=1,ncat ; do j=jsc,jec ; do i=isc,iec
IST%t_surf(i,j,k) = Rad%t_skin(i,j,k)
enddo ; enddo ; enddo
endif
! Determine the sea-ice optical properties.
! These initialization calls for ice-free categories are not really
! needed because these arrays are only used where there is ice.