forked from NOAA-GFDL/SIS2
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathice_model.F90
4284 lines (3798 loc) · 207 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 part of the Flexible Modeling System (FMS). !!
!! !!
!! FMS 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. !!
!! !!
!! FMS 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. !!
!! !!
!! You should have received a copy of the GNU General Public License !!
!! along with FMS; if not, write to: !!
!! Free Software Foundation, Inc. !!
!! 59 Temple Place, Suite 330 !!
!! Boston, MA 02111-1307 USA !!
!! or see: !!
!! http://www.gnu.org/licenses/gpl.txt !!
!! !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
! 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. !
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
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_error_checking, only : chksum, Bchksum, hchksum, uchksum, vchksum
use SIS_error_checking, only : check_redundant_B, check_redundant_C
use SIS_get_input, only : Get_SIS_input, directories
use SIS_sum_output, only : write_ice_statistics, SIS_sum_output_init
use SIS_sum_output, only : accumulate_bottom_input, accumulate_input_1, accumulate_input_2
use MOM_domains, only : pass_var, pass_vector, AGRID, BGRID_NE, CGRID_NE
use MOM_domains, only : fill_symmetric_edges, clone_MOM_domain
use MOM_error_handler, only : SIS_error=>MOM_error, FATAL, WARNING, SIS_mesg=>MOM_mesg
use MOM_file_parser, only : get_param, log_param, log_version, param_file_type
use MOM_file_parser, only : open_param_file, close_param_file
use MOM_string_functions, only : uppercase
use MOM_EOS, only : EOS_type, calculate_density_derivs
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 : 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_LOOP, CLOCK_ROUTINE
use time_manager_mod, only : time_type, time_type_to_real, get_date, get_time
use time_manager_mod, only : set_date, set_time, operator(+), operator(-)
use MOM_time_manager, only : operator(>), operator(*), operator(/), operator(/=)
use 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 constants_mod, only: hlv, hlf, T_0degC=>Tfreeze, grav, STEFAN
use data_override_mod, only : data_override
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, ice_state_type
use ice_type_mod, only : ice_model_restart, dealloc_ice_arrays, dealloc_IST_arrays
use ice_type_mod, only : ice_data_type_register_restarts, ice_state_register_restarts
use ice_type_mod, only : ice_diagnostics_init, ice_stock_pe, check_ice_model_nml
use ice_type_mod, only : ocean_ice_boundary_type, atmos_ice_boundary_type, land_ice_boundary_type
use ice_type_mod, only : ocn_ice_bnd_type_chksum, atm_ice_bnd_type_chksum
use ice_type_mod, only : lnd_ice_bnd_type_chksum, ice_data_type_chksum
use ice_type_mod, only : IST_chksum, Ice_public_type_chksum
use ice_type_mod, only : IST_bounds_check, Ice_public_type_bounds_check
use ice_utils_mod, only : get_avg, post_avg, ice_line, ice_grid_chksum
use SIS_hor_grid_mod, only : SIS_hor_grid_type, set_hor_grid, SIS_hor_grid_end
use SIS_grid_initialize, only : initialize_fixed_SIS_grid
use ice_grid_mod, 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 ice_thm_mod, only: slab_ice_optics, ice_thm_param, ice5lay_temp, ice5lay_resize
use ice_thm_mod, only: MU_TS, TFI, CI, e_to_melt, get_thermo_coefs
use SIS2_ice_thm, only: ice_temp_SIS2, ice_optics_SIS2, SIS2_ice_thm_init, SIS2_ice_thm_end
use SIS2_ice_thm, only: get_SIS2_thermo_coefs, enthalpy_liquid_freeze
use SIS2_ice_thm, only: ice_resize_SIS2, add_frazil_SIS2, rebalance_ice_layers
use SIS2_ice_thm, only: enthalpy_from_TS, enth_from_TS, Temp_from_En_S, Temp_from_Enth_S
use SIS2_ice_thm, only: T_freeze, calculate_T_freeze, enthalpy_liquid, e_to_melt_TS
use ice_dyn_bgrid, only: ice_B_dynamics, ice_B_dyn_init, ice_B_dyn_register_restarts, ice_B_dyn_end
use ice_dyn_cgrid, only: ice_C_dynamics, ice_C_dyn_init, ice_C_dyn_register_restarts, ice_C_dyn_end
use ice_transport_mod, only : ice_transport, ice_transport_init, ice_transport_end
use ice_transport_mod, only : adjust_ice_categories
use ice_bergs, only: icebergs_run, icebergs_init, icebergs_end, icebergs_incr_mass
! ### Eliminate cell_area once flux_exchange.F90 is fixed.
use constants_mod, only : radius, pi
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, ice_model_end, update_ice_model_fast, ice_stock_pe
public :: update_ice_model_slow_up, update_ice_model_slow_dn
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
! ### Eliminate cell_area once flux_exchange.F90 is fixed.
public :: cell_area
integer :: iceClock, iceClock1, iceCLock2, iceCLock3, iceClock4, iceClock5, &
iceClock6, iceClock7, iceClock8, iceClock9, iceClocka, iceClockb, iceClockc
! This is still here as an artefact of an older public interface and should go.
! ### Eliminate cell_area once flux_exchange.F90 is fixed.
real, allocatable, dimension(:,:) :: cell_area ! grid cell area; sphere frac.
contains
!-----------------------------------------------------------------------
!
! Coupler interface to do slow ice processes: dynamics, transport, mass
!
subroutine update_ice_model_slow_dn ( Atmos_boundary, Land_boundary, Ice )
type(atmos_ice_boundary_type), intent(inout) :: Atmos_boundary
type(land_ice_boundary_type), intent(inout) :: Land_boundary
type(ice_data_type), intent(inout) :: Ice
call mpp_clock_begin(iceClock)
call mpp_clock_begin(iceClock2)
call update_ice_model_slow(Ice, Ice%Ice_state, Ice%G, Ice%IG, &
Land_boundary%runoff, Land_boundary%calving, &
Land_boundary%runoff_hflx, Land_boundary%calving_hflx )
call mpp_clock_end(iceClock2)
call mpp_clock_end(iceClock)
end subroutine update_ice_model_slow_dn
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
! sum_top_quantities - sum fluxes for later use by ice/ocean slow physics. !
! Nothing here will be exposed to other modules. !
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
subroutine sum_top_quantities ( Ice, IST, Atmos_boundary_fluxes, flux_u, flux_v, flux_t, flux_q, &
flux_sw_nir_dir, flux_sw_nir_dif, flux_sw_vis_dir, flux_sw_vis_dif,&
flux_lw, lprec, fprec, flux_lh, G, IG)
type(ice_data_type), intent(inout) :: Ice
type(ice_state_type), intent(inout) :: IST
type(coupler_3d_bc_type), intent(inout) :: Atmos_boundary_fluxes
type(SIS_hor_grid_type), intent(inout) :: G
type(ice_grid_type), intent(inout) :: IG
real, dimension(G%isc:G%iec,G%jsc:G%jec,0:IG%CatIce), intent(in) :: &
flux_u, flux_v, flux_t, flux_q, flux_lw, lprec, fprec, flux_lh
real, dimension(G%isc:G%iec,G%jsc:G%jec,0:IG%CatIce), intent(in) :: &
flux_sw_nir_dir, flux_sw_nir_dif, flux_sw_vis_dir, flux_sw_vis_dif
real,dimension(G%isc:G%iec,G%jsc:G%jec) :: tmp
integer :: i, j, k, m, n, i2, j2, k2, isc, iec, jsc, jec, i_off, j_off, ncat
integer :: ind, max_num_fields, next_index
isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; ncat = IG%CatIce
i_off = LBOUND(Ice%albedo_vis_dir,1) - G%isc
j_off = LBOUND(Ice%albedo_vis_dir,2) - G%jsc
if (IST%num_tr_fluxes < 0) then
! Determine how many atmospheric boundary fluxes have been passed in, and
! set up both an indexing array for these and a space to take their average.
! This code is only exercised the first time that sum_top_quantities is called.
IST%num_tr_fluxes = 0
if (Atmos_boundary_fluxes%num_bcs > 0) then
max_num_fields = 0
do n=1,Atmos_boundary_fluxes%num_bcs
IST%num_tr_fluxes = IST%num_tr_fluxes + Atmos_boundary_fluxes%bc(n)%num_fields
max_num_fields = max(max_num_fields, Atmos_boundary_fluxes%bc(n)%num_fields)
enddo
if (IST%num_tr_fluxes > 0) then
allocate(IST%tr_flux_top(SZI_(G), SZJ_(G), 0:IG%CatIce, IST%num_tr_fluxes))
allocate(IST%tr_flux_ocn_top(SZI_(G), SZJ_(G), IST%num_tr_fluxes))
IST%tr_flux_top(:,:,:,:) = 0.0 ; IST%tr_flux_ocn_top(:,:,:) = 0.0
allocate(IST%tr_flux_index(max_num_fields, Atmos_boundary_fluxes%num_bcs))
IST%tr_flux_index(:,:) = -1 ; next_index = 1
do n=1,Atmos_boundary_fluxes%num_bcs ; do m=1,Atmos_boundary_fluxes%bc(n)%num_fields
IST%tr_flux_index(m, n) = next_index ; next_index = next_index + 1
enddo ; enddo
endif
endif
endif
if (IST%avg_count == 0) then
! zero_top_quantities - zero fluxes to begin summing in ice fast physics.
IST%flux_u_top(:,:,:) = 0.0 ; IST%flux_v_top(:,:,:) = 0.0
IST%lwdn(:,:) = 0.0 ; IST%swdn(:,:) = 0.0
IST%flux_t_top(:,:,:) = 0.0 ; IST%flux_q_top(:,:,:) = 0.0
IST%flux_lw_top(:,:,:) = 0.0 ; IST%flux_lh_top(:,:,:) = 0.0
IST%flux_sw_nir_dir_top(:,:,:) = 0.0 ; IST%flux_sw_nir_dif_top(:,:,:) = 0.0
IST%flux_sw_vis_dir_top(:,:,:) = 0.0 ; IST%flux_sw_vis_dif_top(:,:,:) = 0.0
IST%lprec_top(:,:,:) = 0.0 ; IST%fprec_top(:,:,:) = 0.0
if (IST%num_tr_fluxes > 0) IST%tr_flux_top(:,:,:,:) = 0.0
endif
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,ncat,IST,flux_u,flux_v,flux_t, &
!$OMP flux_q,flux_sw_nir_dir,flux_sw_nir_dif, &
!$OMP flux_sw_vis_dir,flux_sw_vis_dif,flux_lw, &
!$OMP lprec,fprec,flux_lh)
do j=jsc,jec ; do k=0,ncat ; do i=isc,iec
IST%flux_u_top(i,j,k) = IST%flux_u_top(i,j,k) + flux_u(i,j,k)
IST%flux_v_top(i,j,k) = IST%flux_v_top(i,j,k) + flux_v(i,j,k)
IST%flux_t_top(i,j,k) = IST%flux_t_top(i,j,k) + flux_t(i,j,k)
IST%flux_q_top(i,j,k) = IST%flux_q_top(i,j,k) + flux_q(i,j,k)
IST%flux_sw_nir_dir_top(i,j,k) = IST%flux_sw_nir_dir_top(i,j,k) + flux_sw_nir_dir(i,j,k)
IST%flux_sw_nir_dif_top(i,j,k) = IST%flux_sw_nir_dif_top(i,j,k) + flux_sw_nir_dif(i,j,k)
IST%flux_sw_vis_dir_top(i,j,k) = IST%flux_sw_vis_dir_top(i,j,k) + flux_sw_vis_dir(i,j,k)
IST%flux_sw_vis_dif_top(i,j,k) = IST%flux_sw_vis_dif_top(i,j,k) + flux_sw_vis_dif(i,j,k)
IST%flux_lw_top(i,j,k) = IST%flux_lw_top(i,j,k) + flux_lw(i,j,k)
IST%lprec_top(i,j,k) = IST%lprec_top(i,j,k) + lprec(i,j,k)
IST%fprec_top(i,j,k) = IST%fprec_top(i,j,k) + fprec(i,j,k)
IST%flux_lh_top(i,j,k) = IST%flux_lh_top(i,j,k) + flux_lh(i,j,k)
enddo ; enddo ; enddo
do n=1,Atmos_boundary_fluxes%num_bcs ; do m=1,Atmos_boundary_fluxes%bc(n)%num_fields
ind = IST%tr_flux_index(m,n)
if (ind < 1) call SIS_error(FATAL, "Bad boundary flux index in sum_top_quantities.")
do k=0,ncat ; do j=jsc,jec ; do i=isc,iec
i2 = i+i_off ; j2 = j+j_off ; k2 = k+1
IST%tr_flux_top(i,j,k,ind) = IST%tr_flux_top(i,j,k,ind) + &
Atmos_boundary_fluxes%bc(n)%field(m)%values(i2,j2,k2)
enddo ; enddo ; enddo
enddo ; enddo
if (IST%id_lwdn > 0) then
call get_avg(flux_lw(:,:,:) + STEFAN*IST%t_surf(isc:iec,jsc:jec,:)**4, &
IST%part_size(isc:iec,jsc:jec,:), tmp(:,:))
do j=jsc,jec ; do i=isc,iec
if (G%mask2dT(i,j)>0.5) IST%lwdn(i,j) = IST%lwdn(i,j) + tmp(i,j)
enddo ; enddo
endif
if (IST%id_swdn > 0) then
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,ncat,G,IST,Ice,i_off,j_off, &
!$OMP flux_sw_vis_dir,flux_sw_vis_dif, &
!$OMP flux_sw_nir_dir,flux_sw_nir_dif) &
!$OMP private(i2,j2,k2)
do j=jsc,jec ; do k=0,ncat ; do i=isc,iec ; if (G%mask2dT(i,j)>0.5) then
i2 = i+i_off ; j2 = j+j_off ; k2 = k+1
IST%swdn(i,j) = IST%swdn(i,j) + IST%part_size(i,j,k) * ( &
(flux_sw_vis_dir(i,j,k)/(1-Ice%albedo_vis_dir(i2,j2,k2)) + &
flux_sw_vis_dif(i,j,k)/(1-Ice%albedo_vis_dif(i2,j2,k2))) + &
(flux_sw_nir_dir(i,j,k)/(1-Ice%albedo_nir_dir(i2,j2,k2)) + &
flux_sw_nir_dif(i,j,k)/(1-Ice%albedo_nir_dif(i2,j2,k2))) )
endif ; enddo ; enddo ; enddo
endif
IST%avg_count = IST%avg_count + 1
end subroutine sum_top_quantities
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
! avg_top_quantities - time average fluxes for ice and ocean slow physics !
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
subroutine avg_top_quantities(Ice, IST, G, IG)
type(ice_data_type), intent(inout) :: Ice
type(ice_state_type), intent(inout) :: IST
type(SIS_hor_grid_type), intent(inout) :: G
type(ice_grid_type), intent(inout) :: IG
real :: u, v, divid, sign
real, dimension(G%isd:G%ied,G%jsd:G%jed) :: tmp2d
integer :: i, j, k, m, n, isc, iec, jsc, jec, ncat
logical :: sent
isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; ncat = IG%CatIce
!
! compute average fluxes
!
if (IST%avg_count == 0) call SIS_error(FATAL,'avg_top_quantities: '//&
'no ocean model fluxes have been averaged')
! Rotate the stress from lat/lon to ocean coordinates and possibly change the
! sign to positive for downward fluxes of positive momentum.
sign = 1.0 ; if (IST%atmos_winds) sign = -1.0
divid = 1.0/real(IST%avg_count)
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,ncat,IST,sign,divid,G) private(u,v)
do j=jsc,jec
do k=0,ncat ; do i=isc,iec
u = IST%flux_u_top(i,j,k) * (sign*divid)
v = IST%flux_v_top(i,j,k) * (sign*divid)
IST%flux_u_top(i,j,k) = u*G%cos_rot(i,j)-v*G%sin_rot(i,j) ! rotate stress from lat/lon
IST%flux_v_top(i,j,k) = v*G%cos_rot(i,j)+u*G%sin_rot(i,j) ! to ocean coordinates
IST%flux_t_top(i,j,k) = IST%flux_t_top(i,j,k) * divid
IST%flux_q_top(i,j,k) = IST%flux_q_top(i,j,k) * divid
IST%flux_sw_nir_dir_top(i,j,k) = IST%flux_sw_nir_dir_top(i,j,k) * divid
IST%flux_sw_nir_dif_top(i,j,k) = IST%flux_sw_nir_dif_top(i,j,k) * divid
IST%flux_sw_vis_dir_top(i,j,k) = IST%flux_sw_vis_dir_top(i,j,k) * divid
IST%flux_sw_vis_dif_top(i,j,k) = IST%flux_sw_vis_dif_top(i,j,k) * divid
IST%flux_lw_top(i,j,k) = IST%flux_lw_top(i,j,k) * divid
IST%fprec_top(i,j,k) = IST%fprec_top(i,j,k) * divid
IST%lprec_top(i,j,k) = IST%lprec_top(i,j,k) * divid
IST%flux_lh_top(i,j,k) = IST%flux_lh_top(i,j,k) * divid
! Convert frost forming atop sea ice into frozen precip.
if ((k>0) .and. (IST%flux_q_top(i,j,k) < 0.0)) then
IST%fprec_top(i,j,k) = IST%fprec_top(i,j,k) - IST%flux_q_top(i,j,k)
IST%flux_q_top(i,j,k) = 0.0
endif
do n=1,IST%num_tr_fluxes
IST%tr_flux_top(i,j,k,n) = IST%tr_flux_top(i,j,k,n) * divid
enddo
enddo ; enddo
do i=isc,iec
IST%lwdn(i,j) = IST%lwdn(i,j)* divid
IST%swdn(i,j) = IST%swdn(i,j)* divid
enddo
enddo
call pass_vector(IST%flux_u_top, IST%flux_v_top, G%Domain, stagger=AGRID)
! Flux diagnostics
!
if (IST%id_sh>0) call post_avg(IST%id_sh, IST%flux_t_top, IST%part_size, &
IST%diag, G=G)
if (IST%id_lh>0) call post_avg(IST%id_lh, IST%flux_lh_top, IST%part_size, &
IST%diag, G=G)
if (IST%id_evap>0) call post_avg(IST%id_evap, IST%flux_q_top, IST%part_size, &
IST%diag, G=G)
if (IST%id_sw>0) then
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,ncat,tmp2d,IST)
do j=jsc,jec
do i=isc,iec ; tmp2d(i,j) = 0.0 ; enddo
do k=0,ncat ; do i=isc,iec
tmp2d(i,j) = tmp2d(i,j) + IST%part_size(i,j,k) * ( &
IST%flux_sw_vis_dir_top(i,j,k) + IST%flux_sw_vis_dif_top(i,j,k) + &
IST%flux_sw_nir_dir_top(i,j,k) + IST%flux_sw_nir_dif_top(i,j,k) )
enddo ; enddo
enddo
call post_data(IST%id_sw, tmp2d, IST%diag)
endif
if (IST%id_lw>0) call post_avg(IST%id_lw, IST%flux_lw_top, &
IST%part_size, IST%diag, G=G)
if (IST%id_snofl>0) call post_avg(IST%id_snofl, IST%fprec_top, &
IST%part_size, IST%diag, G=G)
if (IST%id_rain>0) call post_avg(IST%id_rain, IST%lprec_top, &
IST%part_size, IST%diag, G=G)
if (IST%id_lwdn>0) call post_data(IST%id_lwdn, IST%lwdn, IST%diag)
if (IST%id_swdn>0) call post_data(IST%id_swdn, IST%swdn, IST%diag)
if (IST%id_sw_vis>0) then
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,ncat,tmp2d,IST)
do j=jsc,jec
do i=isc,iec ; tmp2d(i,j) = 0.0 ; enddo
do k=0,ncat ; do i=isc,iec
tmp2d(i,j) = tmp2d(i,j) + IST%part_size(i,j,k) * ( &
IST%flux_sw_vis_dir_top(i,j,k) + IST%flux_sw_vis_dif_top(i,j,k) )
enddo ; enddo
enddo
call post_data(IST%id_sw_vis, tmp2d, IST%diag)
endif
if (IST%id_sw_nir_dir>0) call post_avg(IST%id_sw_nir_dir, IST%flux_sw_nir_dir_top, &
IST%part_size, IST%diag, G=G)
if (IST%id_sw_nir_dif>0) call post_avg(IST%id_sw_nir_dif, IST%flux_sw_nir_dif_top, &
IST%part_size, IST%diag, G=G)
if (IST%id_sw_vis_dir>0) call post_avg(IST%id_sw_vis_dir, IST%flux_sw_vis_dir_top, &
IST%part_size, IST%diag, G=G)
if (IST%id_sw_vis_dif>0) call post_avg(IST%id_sw_vis_dif, IST%flux_sw_vis_dif_top, &
IST%part_size, IST%diag, G=G)
!
! set count to zero and fluxes will be zeroed before the next sum
!
IST%avg_count = 0
end subroutine avg_top_quantities
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
! set_ocean_top_fluxes - Translate 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, G, IG)
type(ice_data_type), intent(inout) :: Ice
type(ice_state_type), intent(inout) :: IST
type(SIS_hor_grid_type), intent(inout) :: G
type(ice_grid_type), intent(inout) :: IG
real :: I_count
integer :: i, j, isc, iec, jsc, jec, m, n, i2, j2, i_off, j_off, ind
isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec
i_off = LBOUND(Ice%flux_t,1) - G%isc ; j_off = LBOUND(Ice%flux_t,2) - G%jsc
if (IST%debug) then
call IST_chksum("Start set_ocean_top_fluxes", IST, G, IG)
call Ice_public_type_chksum("Start set_ocean_top_fluxes", Ice)
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
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,Ice,IST,i_off,j_off) &
!$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) = IST%flux_t_ocn_top(i,j)
Ice%flux_q(i2,j2) = IST%flux_q_ocn_top(i,j)
Ice%flux_sw_vis_dir(i2,j2) = IST%flux_sw_vis_dir_ocn(i,j)
Ice%flux_sw_vis_dif(i2,j2) = IST%flux_sw_vis_dif_ocn(i,j)
Ice%flux_sw_nir_dir(i2,j2) = IST%flux_sw_nir_dir_ocn(i,j)
Ice%flux_sw_nir_dif(i2,j2) = IST%flux_sw_nir_dif_ocn(i,j)
Ice%flux_lw(i2,j2) = IST%flux_lw_ocn_top(i,j)
Ice%flux_lh(i2,j2) = IST%flux_lh_ocn_top(i,j)
Ice%fprec(i2,j2) = IST%fprec_ocn_top(i,j)
Ice%lprec(i2,j2) = IST%lprec_ocn_top(i,j)
enddo ; enddo
if (IST%nudge_sea_ice) 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) + IST%melt_nudge(i,j)
enddo ; enddo
endif
do n=1,Ice%ocean_fluxes%num_bcs ; do m=1,Ice%ocean_fluxes%bc(n)%num_fields
ind = IST%tr_flux_index(m,n)
if (ind < 1) call SIS_error(FATAL, "Bad boundary flux index in set_ocean_top_fluxes.")
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) = IST%tr_flux_ocn_top(i,j,ind)
enddo ; enddo
enddo ; enddo
if (IST%debug) then
call IST_chksum("End set_ocean_top_fluxes", IST, G, IG)
call Ice_public_type_chksum("End set_ocean_top_fluxes", Ice)
endif
end subroutine set_ocean_top_fluxes
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
! finish_ocean_top_stresses - Finish setting the ice-ocean stresses by dividing!
! them through the stresses by the number of times they have been augmented. !
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
subroutine finish_ocean_top_stresses(Ice, IST, G)
type(ice_data_type), intent(inout) :: Ice
type(ice_state_type), intent(inout) :: IST
type(SIS_hor_grid_type), intent(inout) :: G
real :: I_count
if (IST%stress_count > 1) then
I_count = 1.0 / IST%stress_count
Ice%flux_u(:,:) = Ice%flux_u(:,:) * I_count
Ice%flux_v(:,:) = Ice%flux_v(:,:) * I_count
endif
if (IST%debug) then
call Ice_public_type_chksum("finish_ocean_top_stresses", Ice)
endif
end subroutine finish_ocean_top_stresses
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
! set_ocean_top_stress_Bgrid - Calculate the stresses on the ocean integrated !
! across all the thickness categories with the appropriate staggering, and !
! store them in the public ice data type for use by the ocean model. This !
! version of the routine uses wind and ice-ocean stresses on a B-grid. !
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
subroutine set_ocean_top_stress_Bgrid(Ice, IST, windstr_x_water, windstr_y_water, &
str_ice_oce_x, str_ice_oce_y, part_size, G, IG)
type(ice_data_type), intent(inout) :: Ice
type(ice_state_type), intent(inout) :: IST
type(SIS_hor_grid_type), intent(inout) :: G
type(ice_grid_type), intent(inout) :: IG
real, dimension(SZIB_(G),SZJB_(G)), intent(in) :: windstr_x_water, str_ice_oce_x
real, dimension(SZIB_(G),SZJB_(G)), intent(in) :: windstr_y_water, str_ice_oce_y
real, dimension (SZI_(G),SZJ_(G),0:IG%CatIce), intent(in) :: part_size
real :: ps_vel ! part_size interpolated to a velocity point, nondim.
integer :: i, j, k, isc, iec, jsc, jec, ncat, i2, j2, k2, i_off, j_off
isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; ncat = IG%CatIce
i_off = LBOUND(Ice%flux_t,1) - G%isc ; j_off = LBOUND(Ice%flux_t,2) - G%jsc
if (IST%debug) then
call IST_chksum("Start set_ocean_top_stress_Bgrid", IST, G, IG)
call Ice_public_type_chksum("Start set_ocean_top_stress_Bgrid", Ice)
endif
if (IST%stress_count == 0) then
Ice%flux_u(:,:) = 0.0 ; Ice%flux_v(:,:) = 0.0
endif
! Copy and interpolate the ice-ocean stress_Bgrid. This code is slightly
! complicated because there are 3 different staggering options supported.
!$OMP parallel default(none) shared(isc,iec,jsc,jec,ncat,G,Ice,i_off,j_off, &
!$OMP part_size,windstr_x_water,windstr_y_water, &
!$OMP str_ice_oce_x,str_ice_oce_y) &
!$OMP private(i2,j2,ps_vel)
if (Ice%flux_uv_stagger == AGRID) then
!$OMP do
do j=jsc,jec
do i=isc,iec
i2 = i+i_off ; j2 = j+j_off ! Use these to correct for indexing differences.
ps_vel = G%mask2dT(i,j) * part_size(i,j,0)
Ice%flux_u(i2,j2) = Ice%flux_u(i2,j2) + ps_vel * 0.25 * &
((windstr_x_water(I,J) + windstr_x_water(I-1,J-1)) + &
(windstr_x_water(I-1,J) + windstr_x_water(I,J-1)))
Ice%flux_v(i2,j2) = Ice%flux_v(i2,j2) + ps_vel * 0.25 * &
((windstr_y_water(I,J) + windstr_y_water(I-1,J-1)) + &
(windstr_y_water(I-1,J) + windstr_y_water(I,J-1)))
enddo
do k=1,ncat ; do i=isc,iec ; if (G%mask2dT(i,j)>0.5) then
i2 = i+i_off ; j2 = j+j_off ! Use these to correct for indexing differences.
Ice%flux_u(i2,j2) = Ice%flux_u(i2,j2) + part_size(i,j,k) * 0.25 * &
((str_ice_oce_x(I,J) + str_ice_oce_x(I-1,J-1)) + &
(str_ice_oce_x(I-1,J) + str_ice_oce_x(I,J-1)))
Ice%flux_v(i2,j2) = Ice%flux_v(i2,j2) + part_size(i,j,k) * 0.25 * &
((str_ice_oce_y(I,J) + str_ice_oce_y(I-1,J-1)) + &
(str_ice_oce_y(I-1,J) + str_ice_oce_y(I,J-1)))
endif ; enddo ; enddo
enddo
elseif (Ice%flux_uv_stagger == BGRID_NE) then
!$OMP do
do j=jsc,jec
do i=isc,iec
i2 = i+i_off ; j2 = j+j_off ! Use these to correct for indexing differences.
ps_vel = 1.0 ; if (G%mask2dBu(I,J)>0.5) ps_vel = &
0.25*((part_size(i+1,j+1,0) + part_size(i,j,0)) + &
(part_size(i+1,j,0) + part_size(i,j+1,0)) )
Ice%flux_u(i2,j2) = Ice%flux_u(i2,j2) + windstr_x_water(I,J) * ps_vel
Ice%flux_v(i2,j2) = Ice%flux_v(i2,j2) + windstr_y_water(I,J) * ps_vel
enddo
do k=1,ncat ; do i=isc,iec ; if (G%mask2dBu(I,J)>0.5) then
i2 = i+i_off ; j2 = j+j_off ! Use these to correct for indexing differences.
ps_vel = 0.25 * ((part_size(i+1,j+1,k) + part_size(i,j,k)) + &
(part_size(i+1,j,k) + part_size(i,j+1,k)) )
Ice%flux_u(i2,j2) = Ice%flux_u(i2,j2) + str_ice_oce_x(I,J) * ps_vel
Ice%flux_v(i2,j2) = Ice%flux_v(i2,j2) + str_ice_oce_y(I,J) * ps_vel
endif ; enddo ; enddo
enddo
elseif (Ice%flux_uv_stagger == CGRID_NE) then
!$OMP do
do j=jsc,jec
do i=isc,iec
i2 = i+i_off ; j2 = j+j_off ! Use these to correct for indexing differences.
ps_vel = 1.0 ; if (G%mask2dCu(I,j)>0.5) ps_vel = &
0.5*(part_size(i+1,j,0) + part_size(i,j,0))
Ice%flux_u(i2,j2) = Ice%flux_u(i2,j2) + ps_vel * &
0.5 * (windstr_x_water(I,J) + windstr_x_water(I,J-1))
ps_vel = 1.0 ; if (G%mask2dCv(i,J)>0.5) ps_vel = &
0.5*(part_size(i,j+1,0) + part_size(i,j,0))
Ice%flux_v(i2,j2) = Ice%flux_v(i2,j2) + ps_vel * &
0.5 * (windstr_y_water(I,J) + windstr_y_water(I-1,J))
enddo
do k=1,ncat ; do i=isc,iec
i2 = i+i_off ; j2 = j+j_off ! Use these to correct for indexing differences.
if (G%mask2dCu(I,j)>0.5) then
ps_vel = 0.5 * (part_size(i+1,j,k) + part_size(i,j,k))
Ice%flux_u(i2,j2) = Ice%flux_u(i2,j2) + ps_vel * &
0.5 * (str_ice_oce_x(I,J) + str_ice_oce_x(I,J-1))
endif
if (G%mask2dCv(i,J)>0.5) then
ps_vel = 0.5 * (part_size(i,j+1,k) + part_size(i,j,k))
Ice%flux_v(i2,j2) = Ice%flux_v(i2,j2) + ps_vel * &
0.5 * (str_ice_oce_y(I,J) + str_ice_oce_y(I-1,J))
endif
enddo ; enddo
enddo
else
!$OMP single
call SIS_error(FATAL, "set_ocean_top_stress_Bgrid: Unrecognized flux_uv_stagger.")
!$OMP end single
endif
!$OMP end parallel
IST%stress_count = IST%stress_count + 1
if (IST%debug) then
call IST_chksum("End set_ocean_top_stress_Bgrid", IST, G, IG)
call Ice_public_type_chksum("End set_ocean_top_stress_Bgrid", Ice)
endif
end subroutine set_ocean_top_stress_Bgrid
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
! set_ocean_top_stress_Cgrid - Calculate the stresses on the ocean integrated !
! across all the thickness categories with the appropriate staggering, and !
! store them in the public ice data type for use by the ocean model. This !
! version of the routine uses wind and ice-ocean stresses on a C-grid. !
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
subroutine set_ocean_top_stress_Cgrid(Ice, IST, windstr_x_water, windstr_y_water, &
str_ice_oce_x, str_ice_oce_y, part_size, G, IG)
type(ice_data_type), intent(inout) :: Ice
type(ice_state_type), intent(inout) :: IST
type(SIS_hor_grid_type), intent(inout) :: G
type(ice_grid_type), intent(inout) :: IG
real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: windstr_x_water, str_ice_oce_x
real, dimension(SZI_(G),SZJB_(G)), intent(in) :: windstr_y_water, str_ice_oce_y
real, dimension (SZI_(G),SZJ_(G),0:IG%CatIce), intent(in) :: part_size
real :: ps_vel ! part_size interpolated to a velocity point, nondim.
integer :: i, j, k, isc, iec, jsc, jec, ncat, i2, j2, k2, i_off, j_off
isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; ncat = IG%CatIce
i_off = LBOUND(Ice%flux_t,1) - G%isc ; j_off = LBOUND(Ice%flux_t,2) - G%jsc
if (IST%debug) then
call IST_chksum("Start set_ocean_top_stress_Cgrid", IST, G, IG)
call Ice_public_type_chksum("Start set_ocean_top_stress_Cgrid", Ice)
endif
if (IST%stress_count == 0) then
Ice%flux_u(:,:) = 0.0 ; Ice%flux_v(:,:) = 0.0
endif
! Copy and interpolate the ice-ocean stress_Cgrid. This code is slightly
! complicated because there are 3 different staggering options supported.
!$OMP parallel default(none) shared(isc,iec,jsc,jec,ncat,G,Ice,i_off,j_off, &
!$OMP part_size,windstr_x_water,windstr_y_water, &
!$OMP str_ice_oce_x,str_ice_oce_y) &
!$OMP private(i2,j2,ps_vel)
if (Ice%flux_uv_stagger == AGRID) then
!$OMP do
do j=jsc,jec
do i=isc,iec
i2 = i+i_off ; j2 = j+j_off ! Use these to correct for indexing differences.
ps_vel = G%mask2dT(i,j) * part_size(i,j,0)
Ice%flux_u(i2,j2) = Ice%flux_u(i2,j2) + ps_vel * 0.5 * &
(windstr_x_water(I,j) + windstr_x_water(I-1,j))
Ice%flux_v(i2,j2) = Ice%flux_v(i2,j2) + ps_vel * 0.5 * &
(windstr_y_water(I,j) + windstr_y_water(i,J-1))
enddo
do k=1,ncat ; do i=isc,iec ; if (G%mask2dT(i,j)>0.5) then
i2 = i+i_off ; j2 = j+j_off ! Use these to correct for indexing differences.
Ice%flux_u(i2,j2) = Ice%flux_u(i2,j2) + part_size(i,j,k) * 0.5 * &
(str_ice_oce_x(I,j) + str_ice_oce_x(I-1,j))
Ice%flux_v(i2,j2) = Ice%flux_v(i2,j2) + part_size(i,j,k) * 0.5 * &
(str_ice_oce_y(I,j) + str_ice_oce_y(i,J-1))
endif ; enddo ; enddo
enddo
elseif (Ice%flux_uv_stagger == BGRID_NE) then
!$OMP do
do j=jsc,jec
do i=isc,iec
i2 = i+i_off ; j2 = j+j_off ! Use these to correct for indexing differences.
ps_vel = 1.0 ; if (G%mask2dBu(I,J)>0.5) ps_vel = &
0.25*((part_size(i+1,j+1,0) + part_size(i,j,0)) + &
(part_size(i+1,j,0) + part_size(i,j+1,0)) )
! Consider deleting the masks here?
Ice%flux_u(i2,j2) = Ice%flux_u(i2,j2) + ps_vel * G%mask2dBu(I,J) * 0.5 * &
(windstr_x_water(I,j) + windstr_x_water(I,j+1))
Ice%flux_v(i2,j2) = Ice%flux_v(i2,j2) + ps_vel * G%mask2dBu(I,J) * 0.5 * &
(windstr_y_water(I,j) + windstr_y_water(i+1,J))
enddo
do k=1,ncat ; do i=isc,iec ; if (G%mask2dBu(I,J)>0.5) then
i2 = i+i_off ; j2 = j+j_off ! Use these to correct for indexing differences.
ps_vel = 0.25 * ((part_size(i+1,j+1,k) + part_size(i,j,k)) + &
(part_size(i+1,j,k) + part_size(i,j+1,k)) )
Ice%flux_u(i2,j2) = Ice%flux_u(i2,j2) + ps_vel * 0.5 * &
(str_ice_oce_x(I,j) + str_ice_oce_x(I,j+1))
Ice%flux_v(i2,j2) = Ice%flux_v(i2,j2) + ps_vel * 0.5 * &
(str_ice_oce_y(I,j) + str_ice_oce_y(i+1,J))
endif ; enddo ; enddo
enddo
elseif (Ice%flux_uv_stagger == CGRID_NE) then
!$OMP do
do j=jsc,jec
do i=isc,iec
i2 = i+i_off ; j2 = j+j_off ! Use these to correct for indexing differences.
ps_vel = 1.0 ; if (G%mask2dCu(I,j)>0.5) ps_vel = &
0.5*(part_size(i+1,j,0) + part_size(i,j,0))
Ice%flux_u(i2,j2) = Ice%flux_u(i2,j2) + ps_vel * windstr_x_water(I,j)
ps_vel = 1.0 ; if (G%mask2dCv(i,J)>0.5) ps_vel = &
0.5*(part_size(i,j+1,0) + part_size(i,j,0))
Ice%flux_v(i2,j2) = Ice%flux_v(i2,j2) + ps_vel * windstr_y_water(i,J)
enddo
do k=1,ncat ; do i=isc,iec
i2 = i+i_off ; j2 = j+j_off ! Use these to correct for indexing differences.
if (G%mask2dCu(I,j)>0.5) then
ps_vel = 0.5 * (part_size(i+1,j,k) + part_size(i,j,k))
Ice%flux_u(i2,j2) = Ice%flux_u(i2,j2) + ps_vel * str_ice_oce_x(I,j)
endif
if (G%mask2dCv(i,J)>0.5) then
ps_vel = 0.5 * (part_size(i,j+1,k) + part_size(i,j,k))
Ice%flux_v(i2,j2) = Ice%flux_v(i2,j2) + ps_vel * str_ice_oce_y(I,j)
endif
enddo ; enddo
enddo
else
!$OMP single
call SIS_error(FATAL, "set_ocean_top_stress_Cgrid: Unrecognized flux_uv_stagger.")
!$OMP end single
endif
!$OMP end parallel
IST%stress_count = IST%stress_count + 1
if (IST%debug) then
call IST_chksum("End set_ocean_top_stress_Cgrid", IST, G, IG)
call Ice_public_type_chksum("End set_ocean_top_stress_Cgrid", Ice)
endif
end subroutine set_ocean_top_stress_Cgrid
!
! Coupler interface to provide ocean surface data to atmosphere.
!
subroutine update_ice_model_slow_up ( Ocean_boundary, Ice )
type(ocean_ice_boundary_type), intent(inout) :: Ocean_boundary
type(ice_data_type), intent(inout) :: Ice
call mpp_clock_begin(iceClock)
call mpp_clock_begin(iceClock1)
call set_ice_surface_state(Ice, Ice%Ice_state, Ocean_boundary%t, Ocean_boundary%u, Ocean_boundary%v, &
Ocean_boundary%frazil, Ocean_boundary, Ice%G, Ice%IG, &
Ocean_boundary%s, Ocean_boundary%sea_level )
call mpp_clock_end(iceClock1)
call mpp_clock_end(iceClock)
end subroutine update_ice_model_slow_up
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
! set_ice_surface_state - prepare surface state for atmosphere fast physics !
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
subroutine set_ice_surface_state(Ice, IST, t_surf_ice_bot, u_surf_ice_bot, v_surf_ice_bot, &
frazil_ice_bot, OIB, G, IG, s_surf_ice_bot, sea_lev_ice_bot )
type(ice_data_type), intent(inout) :: Ice
type(ice_state_type), intent(inout) :: IST
type(SIS_hor_grid_type), intent(inout) :: G
type(ice_grid_type), intent(inout) :: IG
real, dimension(G%isc:G%iec,G%jsc:G%jec), intent(in) :: t_surf_ice_bot, u_surf_ice_bot
real, dimension(G%isc:G%iec,G%jsc:G%jec), intent(in) :: v_surf_ice_bot, frazil_ice_bot
type(ocean_ice_boundary_type), intent(inout) :: OIB
real, dimension(G%isc:G%iec,G%jsc:G%jec), intent(in) :: s_surf_ice_bot, sea_lev_ice_bot
real, dimension(G%isc:G%iec,G%jsc:G%jec) :: m_ice_tot
real, dimension(G%isc:G%iec,G%jsc:G%jec) :: h_ice_input
real, dimension(G%isc:G%iec,G%jsc:G%jec) :: icec, icec_obs
real, dimension(SZI_(G),SZJ_(G)) :: u_nonsym, v_nonsym
real, dimension(IG%NkIce) :: sw_abs_lay
real :: u, v
real :: area_pt
real :: I_Nk
real :: kg_H_Nk ! The conversion factor from units of H to kg/m2 over Nk.
integer :: i, j, k, m, n, i2, j2, k2, isc, iec, jsc, jec, ncat, i_off, j_off
logical :: sent
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 :: dt_slow ! The ice thermodynamics time step
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
I_Nk = 1.0 / IG%NkIce ; kg_H_Nk = IG%H_to_kg_m2 * I_Nk
dt_slow = time_type_to_real(IST%Time_step_slow)
H_to_m_snow = IG%H_to_kg_m2 / IST%Rho_snow ; H_to_m_ice = IG%H_to_kg_m2 / IST%Rho_ice
! pass ocean state through ice on first partition
if (.not. IST%specified_ice) then ! otherwise, already set by update_ice_model_slow
IST%t_surf(isc:iec,jsc:jec,0) = t_surf_ice_bot(isc:iec,jsc:jec)
endif
if (IST%do_init) then
call get_sea_surface(IST%Time, IST%t_surf(isc:iec,jsc:jec,0), IST%part_size(isc:iec,jsc:jec,0:1), &
h_ice_input )
do j=jsc,jec ; do i=isc,iec
IST%mH_ice(i,j,1) = h_ice_input(i,j)*(IST%Rho_ice*IG%kg_m2_to_H)
enddo ; enddo
! Transfer ice to the correct thickness category. If do_ridging=.false.,
! the first call to ice_redistribute has the same result. At present, all
! tracers are initialized to their default values, and snow is set to 0,
! and so do not need to be updated here.
if (IST%do_ridging) then
do j=jsc,jec ; do i=isc,iec ; if (IST%mH_ice(i,j,1) > IG%mH_cat_bound(1)) then
do k=ncat,2,-1 ; if (IST%mH_ice(i,j,1) > IG%mH_cat_bound(k-1)) then
IST%part_size(i,j,k) = IST%part_size(i,j,1)
IST%part_size(i,j,1) = 0.0
IST%mH_ice(i,j,k) = IST%mH_ice(i,j,1) ; IST%mH_ice(i,j,1) = 0.0
! IST%mH_snow(i,j,k) = IST%mH_snow(i,j,1) ; IST%mH_snow(i,j,1) = 0.0
exit ! from k-loop
endif ; enddo
endif ; enddo ; enddo
endif
call pass_var(IST%part_size, G%Domain, complete=.true. )
call pass_var(IST%mH_ice, G%Domain, complete=.true. )
IST%do_init = .false.
endif
! Any special first-time initialization must be completed before this point.
IST%first_time = .false.
if (IST%bounds_check) &
call IST_bounds_check(IST, G, IG, "Start of set_ice_surface_state")
if (IST%debug) then
call IST_chksum("Start set_ice_surface_state", IST, G, IG)
call Ice_public_type_chksum("Start set_ice_surface_state", Ice)
call chksum(u_surf_ice_bot(isc:iec,jsc:jec), "Start IB2IT u_surf_ice_bot")
call chksum(v_surf_ice_bot(isc:iec,jsc:jec), "Start IB2IT v_surf_ice_bot")
endif
! Transfer the ocean state for extra tracer fluxes.
do n=1,OIB%fields%num_bcs ; do m=1,OIB%fields%bc(n)%num_fields
Ice%ocean_fields%bc(n)%field(m)%values(:,:,1) = OIB%fields%bc(n)%field(m)%values
enddo ; enddo
m_ice_tot(:,:) = 0.0
! ###Eliminate ice_mask once flux_exchange.F90 is fixed.
Ice%ice_mask(:,:,1) = .false.
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,G,IST,t_surf_ice_bot, &
!$OMP s_surf_ice_bot,frazil_ice_bot,sea_lev_ice_bot, &
!$OMP ncat,m_ice_tot,Ice,i_off,j_off) &
!$OMP private(i2,j2,k2)
do j=jsc,jec
do i=isc,iec
IST%t_ocn(i,j) = t_surf_ice_bot(i,j) - T_0degC
IST%s_surf(i,j) = s_surf_ice_bot(i,j)
IST%frazil(i,j) = frazil_ice_bot(i,j)
IST%sea_lev(i,j) = sea_lev_ice_bot(i,j)
enddo
do k=1,ncat ; do i=isc,iec
IST%tmelt(i,j,k) = 0.0 ; IST%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
IST%bheat(i,j) = IST%kmelt*(IST%t_ocn(i,j) - T_Freeze(IST%s_surf(i,j), IST%ITV))
else
IST%bheat(i,j) = 0.0
endif
enddo
! ###Eliminate ice_mask once flux_exchange.F90 is fixed.
do k=1,ncat ; do i=isc,iec
i2 = i+i_off ; j2 = j+j_off ; k2 = k+1
Ice%ice_mask(i2,j2,k2) = (IST%mH_ice(i,j,k) > 0.0)
enddo ; enddo
enddo
if (IST%slab_ice) then
IST%sw_abs_sfc(:,:,:) = 0.0 ; IST%sw_abs_snow(:,:,:) = 0.0
IST%sw_abs_ice(:,:,:,:) = 0.0 ; IST%sw_abs_ocn(:,:,:) = 0.0
IST%sw_abs_int(:,:,:) = 0.0
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,ncat,IST,Ice,i_off,j_off, &
!$OMP H_to_m_snow,H_to_m_ice) &
!$OMP private(i2,j2,k2)
do j=jsc,jec ; do k=1,ncat ; do i=isc,iec ; if (IST%mH_ice(i,j,k) > 0.0) then
i2 = i+i_off ; j2 = j+j_off ; k2 = k+1
call slab_ice_optics(IST%mH_snow(i,j,k)*H_to_m_snow, IST%mH_ice(i,j,k)*H_to_m_ice, &
IST%t_surf(i,j,k)-T_0degC, T_Freeze(IST%s_surf(i,j),IST%ITV), &
Ice%albedo(i2,j2,k2))
Ice%albedo_vis_dir(i2,j2,k2) = Ice%albedo(i2,j2,k2)
Ice%albedo_vis_dif(i2,j2,k2) = Ice%albedo(i2,j2,k2)
Ice%albedo_nir_dir(i2,j2,k2) = Ice%albedo(i2,j2,k2)
Ice%albedo_nir_dif(i2,j2,k2) = Ice%albedo(i2,j2,k2)
endif ; enddo ; enddo ; enddo
else
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,ncat,IST,Ice,G,IG,i_off,j_off, &
!$OMP H_to_m_snow,H_to_m_ice) &
!$OMP private(i2,j2,k2,sw_abs_lay)
do j=jsc,jec ; do k=1,ncat ; do i=isc,iec ; if (IST%mH_ice(i,j,k) > 0.0) then
i2 = i+i_off ; j2 = j+j_off ; k2 = k+1
call ice_optics_SIS2(IST%mH_snow(i,j,k)*H_to_m_snow, IST%mH_ice(i,j,k)*H_to_m_ice, &
IST%t_surf(i,j,k)-T_0degC, T_Freeze(IST%s_surf(i,j),IST%ITV), IG%NkIce, &
Ice%albedo_vis_dir(i2,j2,k2), Ice%albedo_vis_dif(i2,j2,k2), &
Ice%albedo_nir_dir(i2,j2,k2), Ice%albedo_nir_dif(i2,j2,k2), &
IST%sw_abs_sfc(i,j,k), IST%sw_abs_snow(i,j,k), &
sw_abs_lay, IST%sw_abs_ocn(i,j,k), IST%sw_abs_int(i,j,k), &
IST%ice_thm_CSp, coszen_in=IST%coszen(i,j))
do m=1,IG%NkIce ; IST%sw_abs_ice(i,j,k,m) = sw_abs_lay(m) ; enddo
!Niki: Is the following correct for diagnostics?
! Probably this calculation of the "average" albedo should be replaced
! with a calculation that weights the averaging by the fraction of the
! shortwave radiation in each wavelength and orientation band. However,
! since this is only used for diagnostic purposes, making this change
! might not be too urgent. -RWH
Ice%albedo(i2,j2,k2)=(Ice%albedo_vis_dir(i2,j2,k2)+Ice%albedo_nir_dir(i2,j2,k2)&
+Ice%albedo_vis_dif(i2,j2,k2)+Ice%albedo_nir_dif(i2,j2,k2))/4
endif ; enddo ; enddo ; enddo
endif
if (Ice%flux_uv_stagger == AGRID) then
u_nonsym(:,:) = 0.0 ; v_nonsym(:,:) = 0.0
do j=jsc,jec ; do i=isc,iec
u_nonsym(i,j) = u_surf_ice_bot(i,j) ; v_nonsym(i,j) = v_surf_ice_bot(i,j)
enddo ; enddo
call pass_vector(u_nonsym, v_nonsym, G%Domain_aux, stagger=AGRID)
if (associated(IST%u_ocn) .and. associated(IST%v_ocn)) then
do J=jsc-1,jec ; do I=isc-1,iec
IST%u_ocn(I,J) = 0.25*((u_nonsym(i,j) + u_nonsym(i+1,j+1)) + &
(u_nonsym(i+1,j) + u_nonsym(i,j+1)))
IST%v_ocn(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
call pass_vector(IST%u_ocn, IST%v_ocn, G%Domain, stagger=BGRID_NE)
endif
if (associated(IST%u_ocn_C) .and. associated(IST%v_ocn_C)) then
do j=jsc,jec ; do I=isc-1,iec
IST%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
IST%v_ocn_C(i,J) = 0.5*(v_nonsym(i,j) + v_nonsym(i,j+1))
enddo ; enddo
call pass_vector(IST%u_ocn_C, IST%v_ocn_C, G%Domain, stagger=CGRID_NE)
endif
elseif (OIB%stagger == BGRID_NE) then
if (IST%Cgrid_dyn) then
u_nonsym(:,:) = 0.0 ; v_nonsym(:,:) = 0.0
do j=jsc,jec ; do i=isc,iec
u_nonsym(i,j) = u_surf_ice_bot(i,j) ; v_nonsym(i,j) = v_surf_ice_bot(i,j)
enddo ; enddo
call pass_vector(u_nonsym, v_nonsym, G%Domain_aux, stagger=BGRID_NE)
do j=jsc,jec ; do I=isc-1,iec
IST%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
IST%v_ocn_C(i,J) = 0.5*(v_nonsym(I,J) + v_nonsym(I-1,J))
enddo ; enddo
call pass_vector(IST%u_ocn_C, IST%v_ocn_C, G%Domain, stagger=CGRID_NE)
else
do J=jsc,jec ; do I=isc,iec
IST%u_ocn(I,J) = u_surf_ice_bot(I,J) ! need under-ice current
IST%v_ocn(I,J) = v_surf_ice_bot(I,J) ! for water drag term
enddo ; enddo
if (G%symmetric) &
call fill_symmetric_edges(IST%u_ocn, IST%v_ocn, G%Domain, stagger=BGRID_NE)
call pass_vector(IST%u_ocn, IST%v_ocn, G%Domain, stagger=BGRID_NE)
endif
elseif (OIB%stagger == CGRID_NE) then
if (IST%Cgrid_dyn) then
do j=jsc,jec ; do I=isc,iec
IST%u_ocn_C(I,j) = u_surf_ice_bot(I,j)
enddo ; enddo
do J=jsc,jec ; do i=isc,iec
IST%v_ocn_C(i,J) = v_surf_ice_bot(I,j)
enddo ; enddo
if (G%symmetric) &
call fill_symmetric_edges(IST%u_ocn_C, IST%v_ocn_C, G%Domain, stagger=CGRID_NE)
call pass_vector(IST%u_ocn_C, IST%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
u_nonsym(I,j) = u_surf_ice_bot(I,j) ; v_nonsym(i,J) = v_surf_ice_bot(i,J)
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
IST%u_ocn(I,J) = 0.5*(u_nonsym(I,j) + u_nonsym(I,j+1))
IST%v_ocn(I,J) = 0.5*(v_nonsym(i,J) + v_nonsym(i+1,J))
enddo ; enddo
call pass_vector(IST%u_ocn, IST%v_ocn, G%Domain, stagger=BGRID_NE)
endif
else
call SIS_error(FATAL, "set_ice_surface_state: Unrecognized OIB%stagger.")
endif
call pass_var(IST%sea_lev, G%Domain)
if (IST%debug) then
if (associated(IST%u_ocn) .and. associated(IST%v_ocn)) then
call chksum(IST%u_ocn(isc:iec,jsc:jec), "Post-pass IST%u_ocn(0,0)")
call chksum(IST%v_ocn(isc:iec,jsc:jec), "Post-pass IST%v_ocn(0,0)")
endif
if (associated(IST%u_ocn_C) .and. associated(IST%v_ocn_C)) then
call chksum(IST%u_ocn_C(isc:iec,jsc:jec), "Post-pass IST%u_ocn_C(0,0)")
call chksum(IST%v_ocn_C(isc:iec,jsc:jec), "Post-pass IST%v_ocn_C(0,0)")
endif
endif
if (IST%bounds_check) &
call IST_bounds_check(IST, G, IG, "Midpoint set_ice_surface_state")
! Copy the surface temperatures into the externally visible data type.