forked from NOAA-GFDL/icebergs
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathicebergs_io.F90
1761 lines (1555 loc) · 73.5 KB
/
icebergs_io.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
!> Handles reading/writing of restart files and trajectory-based diagnostic files
module ice_bergs_io
! This file is part of NOAA-GFDL/icebergs. See LICENSE.md for the license.
use mpp_domains_mod, only: domain2D
use mpp_domains_mod, only: mpp_domain_is_tile_root_pe,mpp_get_domain_tile_root_pe
use mpp_domains_mod, only: mpp_get_tile_pelist,mpp_get_tile_npes,mpp_get_io_domain,mpp_get_tile_id
use mpp_mod, only: mpp_npes, mpp_pe, mpp_root_pe, mpp_sum, mpp_min, mpp_max, NULL_PE
use mpp_mod, only: mpp_send, mpp_recv, mpp_gather, mpp_chksum
use mpp_mod, only: COMM_TAG_11, COMM_TAG_12, COMM_TAG_13, COMM_TAG_14
use fms_mod, only: stdlog, stderr, error_mesg, FATAL, WARNING, NOTE
use fms_mod, only: field_exist, file_exist, read_data, write_data
use fms_io_mod, only: get_instance_filename
use fms_io_mod, only : save_restart, restart_file_type, free_restart_type, set_meta_global
use fms_io_mod, only : register_restart_axis, register_restart_field, set_domain, nullify_domain
use fms_io_mod, only : read_unlimited_axis =>read_compressed, field_exist, get_field_size
use mpp_mod, only : mpp_clock_begin, mpp_clock_end, mpp_clock_id
use mpp_mod, only : CLOCK_COMPONENT, CLOCK_SUBCOMPONENT, CLOCK_LOOP
use fms_mod, only : clock_flag_default
use time_manager_mod, only: time_type, get_date, get_time, set_date, operator(-)
use ice_bergs_framework, only: icebergs_gridded, xyt, iceberg, icebergs, buffer, bond
use ice_bergs_framework, only: pack_traj_into_buffer2,unpack_traj_from_buffer2
use ice_bergs_framework, only: find_cell,find_cell_by_search,count_bergs,is_point_in_cell,pos_within_cell,append_posn
use ice_bergs_framework, only: count_bonds, form_a_bond, find_individual_iceberg
use ice_bergs_framework, only: push_posn
use ice_bergs_framework, only: add_new_berg_to_list,destroy_iceberg
use ice_bergs_framework, only: increase_ibuffer,grd_chksum2,grd_chksum3
use ice_bergs_framework, only: sum_mass,sum_heat,bilin
!params !Niki: write a subroutine to get these
use ice_bergs_framework, only: nclasses, buffer_width, buffer_width_traj
use ice_bergs_framework, only: verbose, really_debug, debug, restart_input_dir,make_calving_reproduce
use ice_bergs_framework, only: ignore_ij_restart, use_slow_find,generate_test_icebergs,print_berg
use ice_bergs_framework, only: force_all_pes_traj
use ice_bergs_framework, only: check_for_duplicates_in_parallel
use ice_bergs_framework, only: split_id, id_from_2_ints, generate_id
implicit none ; private
include 'netcdf.inc'
public ice_bergs_io_init
public read_restart_bergs, write_restart, write_trajectory
public read_restart_calving, read_restart_bonds
public read_ocean_depth
!Local Vars
integer, parameter :: file_format_major_version=0
integer, parameter :: file_format_minor_version=1
!I/O vars
type(domain2d), pointer, save :: io_domain=>NULL()
integer, save :: io_tile_id(1), io_tile_root_pe, io_npes
integer, allocatable,save :: io_tile_pelist(:)
logical :: is_io_tile_root_pe = .true.
integer :: clock_trw,clock_trp
#ifdef _FILE_VERSION
character(len=128) :: version = _FILE_VERSION
#else
character(len=128) :: version = 'unknown'
#endif
contains
!> Initialize parallel i/o
subroutine ice_bergs_io_init(bergs, io_layout)
type(icebergs), pointer :: bergs !< Icebergs container
integer, intent(in) :: io_layout(2) !< Decomposition of i/o processors
integer :: np
integer :: stdlogunit, stderrunit
! Get the stderr and stdlog unit numbers
stderrunit=stderr()
stdlogunit=stdlog()
write(stdlogunit,*) "ice_bergs_framework: "//trim(version)
!I/O layout init
io_tile_id=-1
io_domain => mpp_get_io_domain(bergs%grd%domain)
if(associated(io_domain)) then
io_tile_id = mpp_get_tile_id(io_domain)
is_io_tile_root_pe = mpp_domain_is_tile_root_pe(io_domain)
io_tile_root_pe = mpp_get_domain_tile_root_pe(io_domain)
np=mpp_get_tile_npes(io_domain)
allocate(io_tile_pelist(np))
call mpp_get_tile_pelist(io_domain,io_tile_pelist)
io_npes = io_layout(1)*io_layout(2)
endif
clock_trw=mpp_clock_id( 'Icebergs-traj write', flags=clock_flag_default, grain=CLOCK_SUBCOMPONENT )
clock_trp=mpp_clock_id( 'Icebergs-traj prepare', flags=clock_flag_default, grain=CLOCK_SUBCOMPONENT )
end subroutine ice_bergs_io_init
!> Write an iceberg restart file
subroutine write_restart(bergs, time_stamp)
! Arguments
type(icebergs), pointer :: bergs !< Icebergs container
character(len=*), intent(in), optional :: time_stamp !< Timestamp for restart file
! Local variables
type(bond), pointer :: current_bond
integer :: i,j,id
character(len=35) :: filename
character(len=35) :: filename_bonds
type(iceberg), pointer :: this=>NULL()
integer :: stderrunit
!I/O vars
type(restart_file_type) :: bergs_restart
type(restart_file_type) :: bergs_bond_restart
type(restart_file_type) :: calving_restart
integer :: nbergs, nbonds
integer :: n_static_bergs
logical :: check_bond_quality
type(icebergs_gridded), pointer :: grd
real, allocatable, dimension(:) :: lon, &
lat, &
uvel, &
vvel, &
mass, &
axn, &
ayn, &
bxn, &
byn, &
thickness, &
width, &
length, &
start_lon, &
start_lat, &
start_day, &
start_mass, &
mass_scaling, &
mass_of_bits, &
static_berg, &
heat_density
integer, allocatable, dimension(:) :: ine, &
jne, &
id_cnt, &
id_ij, &
start_year, &
first_id_cnt, &
other_id_cnt, &
first_id_ij, &
other_id_ij, &
first_berg_jne, &
first_berg_ine, &
other_berg_jne, &
other_berg_ine
integer :: grdi, grdj
! Get the stderr unit number
stderrunit=stderr()
! For convenience
grd=>bergs%grd
!First add the bergs on the io_tile_root_pe (if any) to the I/O list
nbergs = 0
do grdj = bergs%grd%jsc,bergs%grd%jec ; do grdi = bergs%grd%isc,bergs%grd%iec
this=>bergs%list(grdi,grdj)%first
do while (associated(this))
nbergs = nbergs +1
this=>this%next
enddo
enddo ; enddo
allocate(lon(nbergs))
allocate(lat(nbergs))
allocate(uvel(nbergs))
allocate(vvel(nbergs))
allocate(mass(nbergs))
allocate(axn(nbergs)) !Alon
allocate(ayn(nbergs)) !Alon
allocate(bxn(nbergs)) !Alon
allocate(byn(nbergs)) !Alon
allocate(thickness(nbergs))
allocate(width(nbergs))
allocate(length(nbergs))
allocate(start_lon(nbergs))
allocate(start_lat(nbergs))
allocate(start_day(nbergs))
allocate(start_mass(nbergs))
allocate(mass_scaling(nbergs))
allocate(mass_of_bits(nbergs))
allocate(heat_density(nbergs))
allocate(static_berg(nbergs))
allocate(ine(nbergs))
allocate(jne(nbergs))
allocate(start_year(nbergs))
allocate(id_cnt(nbergs))
allocate(id_ij(nbergs))
filename = trim("icebergs.res.nc")
call set_domain(bergs%grd%domain)
call register_restart_axis(bergs_restart,filename,'i',nbergs)
call set_meta_global(bergs_restart,'file_format_major_version',ival=(/file_format_major_version/))
call set_meta_global(bergs_restart,'file_format_minor_version',ival=(/file_format_minor_version/))
call set_meta_global(bergs_restart,'time_axis',ival=(/0/))
!Now start writing in the io_tile_root_pe if there are any bergs in the I/O list
! Define Variables
id = register_restart_field(bergs_restart,filename,'lon',lon,longname='longitude',units='degrees_E')
id = register_restart_field(bergs_restart,filename,'lat',lat,longname='latitude',units='degrees_N')
id = register_restart_field(bergs_restart,filename,'uvel',uvel,longname='zonal velocity',units='m/s')
id = register_restart_field(bergs_restart,filename,'vvel',vvel,longname='meridional velocity',units='m/s')
id = register_restart_field(bergs_restart,filename,'mass',mass,longname='mass',units='kg')
if (.not. bergs%Runge_not_Verlet) then
id = register_restart_field(bergs_restart,filename,'axn',axn,longname='explicit zonal acceleration',units='m/s^2')
id = register_restart_field(bergs_restart,filename,'ayn',ayn,longname='explicit meridional acceleration',units='m/s^2')
id = register_restart_field(bergs_restart,filename,'bxn',bxn,longname='inplicit zonal acceleration',units='m/s^2')
id = register_restart_field(bergs_restart,filename,'byn',byn,longname='implicit meridional acceleration',units='m/s^2')
endif
id = register_restart_field(bergs_restart,filename,'ine',ine,longname='i index',units='none')
id = register_restart_field(bergs_restart,filename,'jne',jne,longname='j index',units='none')
id = register_restart_field(bergs_restart,filename,'thickness',thickness,longname='thickness',units='m')
id = register_restart_field(bergs_restart,filename,'width',width,longname='width',units='m')
id = register_restart_field(bergs_restart,filename,'length',length,longname='length',units='m')
id = register_restart_field(bergs_restart,filename,'start_lon',start_lon, &
longname='longitude of calving location',units='degrees_E')
id = register_restart_field(bergs_restart,filename,'start_lat',start_lat, &
longname='latitude of calving location',units='degrees_N')
id = register_restart_field(bergs_restart,filename,'start_year',start_year, &
longname='calendar year of calving event', units='years')
id = register_restart_field(bergs_restart,filename,'id_cnt',id_cnt, &
longname='counter component of iceberg id', units='dimensionless')
id = register_restart_field(bergs_restart,filename,'id_ij',id_ij, &
longname='position component of iceberg id', units='dimensionless')
id = register_restart_field(bergs_restart,filename,'start_day',start_day, &
longname='year day of calving event',units='days')
id = register_restart_field(bergs_restart,filename,'start_mass',start_mass, &
longname='initial mass of calving berg',units='kg')
id = register_restart_field(bergs_restart,filename,'mass_scaling',mass_scaling, &
longname='scaling factor for mass of calving berg',units='none')
id = register_restart_field(bergs_restart,filename,'mass_of_bits',mass_of_bits, &
longname='mass of bergy bits',units='kg')
id = register_restart_field(bergs_restart,filename,'heat_density',heat_density, &
longname='heat density',units='J/kg')
!Checking if any icebergs are static in order to decide whether to save static_berg
n_static_bergs = 0
do grdj = bergs%grd%jsc,bergs%grd%jec ; do grdi = bergs%grd%isc,bergs%grd%iec
this=>bergs%list(grdi,grdj)%first
do while (associated(this))
n_static_bergs=n_static_bergs+this%static_berg
this=>this%next
enddo
enddo ; enddo
call mpp_sum(n_static_bergs)
if (n_static_bergs .gt. 0) &
id = register_restart_field(bergs_restart,filename,'static_berg',static_berg, &
longname='static_berg',units='dimensionless')
! Write variables
i = 0
do grdj = bergs%grd%jsc,bergs%grd%jec ; do grdi = bergs%grd%isc,bergs%grd%iec
this=>bergs%list(grdi,grdj)%first
do while(associated(this))
i = i + 1
lon(i) = this%lon; lat(i) = this%lat
uvel(i) = this%uvel; vvel(i) = this%vvel
ine(i) = this%ine; jne(i) = this%jne
mass(i) = this%mass; thickness(i) = this%thickness
axn(i) = this%axn; ayn(i) = this%ayn !Added by Alon
bxn(i) = this%bxn; byn(i) = this%byn !Added by Alon
width(i) = this%width; length(i) = this%length
start_lon(i) = this%start_lon; start_lat(i) = this%start_lat
start_year(i) = this%start_year; start_day(i) = this%start_day
start_mass(i) = this%start_mass; mass_scaling(i) = this%mass_scaling
static_berg(i) = this%static_berg
call split_id(this%id, id_cnt(i), id_ij(i))
mass_of_bits(i) = this%mass_of_bits; heat_density(i) = this%heat_density
this=>this%next
enddo
enddo ; enddo
call save_restart(bergs_restart, time_stamp)
call free_restart_type(bergs_restart)
deallocate( &
lon, &
lat, &
uvel, &
vvel, &
mass, &
axn, &
ayn, &
bxn, &
byn, &
thickness, &
width, &
length, &
start_lon, &
start_lat, &
start_day, &
start_mass, &
mass_scaling, &
mass_of_bits, &
static_berg, &
heat_density )
deallocate( &
ine, &
jne, &
id_cnt, &
id_ij, &
start_year )
call nullify_domain()
!########## Creating bond restart file ######################
!Allocating restart memory for bond related variables.
nbonds=0
if (bergs%iceberg_bonds_on) then
check_bond_quality=.true.
call count_bonds(bergs, nbonds,check_bond_quality)
allocate(first_id_cnt(nbonds))
allocate(other_id_cnt(nbonds))
allocate(first_id_ij(nbonds))
allocate(other_id_ij(nbonds))
allocate(first_berg_ine(nbonds))
allocate(first_berg_jne(nbonds))
allocate(other_berg_ine(nbonds))
allocate(other_berg_jne(nbonds))
call get_instance_filename("bonds_iceberg.res.nc", filename_bonds)
call set_domain(bergs%grd%domain)
call register_restart_axis(bergs_bond_restart,filename,'i',nbonds)
call set_meta_global(bergs_bond_restart,'file_format_major_version',ival=(/file_format_major_version/))
call set_meta_global(bergs_bond_restart,'file_format_minor_version',ival=(/file_format_minor_version/))
call set_meta_global(bergs_bond_restart,'time_axis',ival=(/0/))
!Now start writing in the io_tile_root_pe if there are any bergs in the I/O list
id = register_restart_field(bergs_bond_restart,filename_bonds,'first_berg_ine',first_berg_ine,longname='iceberg ine of first berg in bond',units='dimensionless')
id = register_restart_field(bergs_bond_restart,filename_bonds,'first_berg_jne',first_berg_jne,longname='iceberg jne of first berg in bond',units='dimensionless')
id = register_restart_field(bergs_bond_restart,filename_bonds,'first_id_cnt',first_id_cnt,longname='counter component of iceberg id first berg in bond',units='dimensionless')
id = register_restart_field(bergs_bond_restart,filename_bonds,'first_id_ij',first_id_ij,longname='position component of iceberg id first berg in bond',units='dimensionless')
id = register_restart_field(bergs_bond_restart,filename_bonds,'other_berg_ine',other_berg_ine,longname='iceberg ine of second berg in bond',units='dimensionless')
id = register_restart_field(bergs_bond_restart,filename_bonds,'other_berg_jne',other_berg_jne,longname='iceberg jne of second berg in bond',units='dimensionless')
id = register_restart_field(bergs_bond_restart,filename_bonds,'other_id_cnt',other_id_cnt,longname='counter component of iceberg id second berg in bond',units='dimensionless')
id = register_restart_field(bergs_bond_restart,filename_bonds,'other_id_ij',other_id_ij,longname='position component of iceberg id second berg in bond',units='dimensionless')
! Write variables
i = 0
do grdj = bergs%grd%jsc,bergs%grd%jec ; do grdi = bergs%grd%isc,bergs%grd%iec
this=>bergs%list(grdi,grdj)%first
do while(associated(this)) !Loops over all bergs
current_bond=>this%first_bond
do while (associated(current_bond)) ! loop over all bonds
i = i + 1
first_berg_ine(i)=this%ine
first_berg_jne(i)=this%jne
call split_id( this%id, first_id_cnt(i), first_id_ij(i) )
call split_id( current_bond%other_id, other_id_cnt(i), other_id_ij(i) )
other_berg_ine(i)=current_bond%other_berg%ine
other_berg_jne(i)=current_bond%other_berg%jne
current_bond=>current_bond%next_bond
enddo !End of loop over bonds
this=>this%next
enddo!End of loop over bergs
enddo; enddo !End of loop over grid
call save_restart(bergs_bond_restart, time_stamp)
call free_restart_type(bergs_bond_restart)
deallocate(first_id_cnt, &
other_id_cnt, &
first_id_ij, &
other_id_ij, &
first_berg_ine, &
first_berg_jne, &
other_berg_ine, &
other_berg_jne )
call nullify_domain()
endif
! Write stored ice
filename='calving.res.nc'
if (verbose.and.mpp_pe().eq.mpp_root_pe()) write(stderrunit,'(2a)') 'diamonds, write_restart: writing ',filename
call grd_chksum3(bergs%grd, bergs%grd%stored_ice, 'write stored_ice')
call grd_chksum2(bergs%grd, bergs%grd%stored_heat, 'write stored_heat')
if (bergs%tau_calving>0.) then
call grd_chksum2(bergs%grd, bergs%grd%rmean_calving, 'write mean calving')
call grd_chksum2(bergs%grd, bergs%grd%rmean_calving_hflx, 'write mean calving_hflx')
endif
call set_domain(bergs%grd%domain)
id = register_restart_field(calving_restart,filename,'stored_ice',bergs%grd%stored_ice,longname='STORED_ICE',units='none')
id = register_restart_field(calving_restart,filename,'stored_heat',bergs%grd%stored_heat,longname='STORED_HEAT',units='none')
id = register_restart_field(calving_restart,filename,'iceberg_counter_grd',bergs%grd%iceberg_counter_grd,longname='ICEBERG_COUNTER_GRD',units='none')
if (bergs%tau_calving>0.) then
id = register_restart_field(calving_restart,filename,'rmean_calving',bergs%grd%rmean_calving,longname='RMEAN_CALVING',units='none')
id = register_restart_field(calving_restart,filename,'rmean_calving_hflx',bergs%grd%rmean_calving_hflx,longname='RMEAN_CALVING_HFLX',units='none')
endif
call save_restart(calving_restart, time_stamp)
call free_restart_type(calving_restart)
end subroutine write_restart
!> Find the last berg in a linked list.
function last_berg(berg)
! Arguments
type(iceberg), pointer :: berg !< Pointer to an iceberg
type(iceberg), pointer :: last_berg
! Local variables
last_berg=>berg
do while (associated(last_berg%next))
last_berg=>last_berg%next
enddo
end function last_berg
!> Read a real value from a file and optionally return a default value if variable is missing
real function get_real_from_file(ncid, varid, k, value_if_not_in_file)
integer, intent(in) :: ncid, varid, k
real, optional :: value_if_not_in_file
if (varid<1.and.present(value_if_not_in_file)) then
get_real_from_file=value_if_not_in_file
else
get_real_from_file=get_double(ncid, ncid, k)
endif
end function get_real_from_file
!> Read an iceberg restart file
subroutine read_restart_bergs(bergs,Time)
! Arguments
type(icebergs), pointer :: bergs !< Icebergs container
type(time_type), intent(in) :: Time !< Model time
! Local variables
integer :: k, siz(4), nbergs_in_file, nbergs_read
logical :: lres, found_restart, found, replace_iceberg_num
logical :: explain
logical :: multiPErestart ! Not needed with new restart read; currently kept for compatibility
real :: lon0, lon1, lat0, lat1
real :: pos_is_good, pos_is_good_all_pe
character(len=33) :: filename, filename_base
type(icebergs_gridded), pointer :: grd
type(iceberg) :: localberg ! NOT a pointer but an actual local variable
integer :: stderrunit, i, j, cnt, ij
real, allocatable, dimension(:) :: lon, &
lat, &
uvel, &
vvel, &
mass, &
axn, &
ayn, &
bxn, &
byn, &
thickness, &
width, &
length, &
start_lon, &
start_lat, &
start_day, &
start_mass, &
mass_scaling, &
mass_of_bits, &
static_berg, &
heat_density
integer, allocatable, dimension(:) :: ine, &
jne, &
iceberg_num,&
id_cnt, &
id_ij, &
start_year
!integer, allocatable, dimension(:,:) :: iceberg_counter_grd
! Get the stderr unit number
stderrunit=stderr()
! For convenience
grd=>bergs%grd
! Zero out nbergs_in_file
nbergs_in_file = 0
filename_base=trim(restart_input_dir)//'icebergs.res.nc'
found_restart = find_restart_file(filename_base, filename, multiPErestart, io_tile_id(1))
if (found_restart) then
filename = filename_base
call get_field_size(filename,'i',siz, field_found=found, domain=bergs%grd%domain)
nbergs_in_file = siz(1)
replace_iceberg_num = field_exist(filename, 'iceberg_num') ! True if using 32-bit iceberg_num in restart file
allocate(lon(nbergs_in_file))
allocate(lat(nbergs_in_file))
allocate(uvel(nbergs_in_file))
allocate(vvel(nbergs_in_file))
allocate(mass(nbergs_in_file))
allocate(axn(nbergs_in_file)) !Alon
allocate(ayn(nbergs_in_file)) !Alon
allocate(bxn(nbergs_in_file)) !Alon
allocate(byn(nbergs_in_file)) !Alon
allocate(thickness(nbergs_in_file))
allocate(width(nbergs_in_file))
allocate(length(nbergs_in_file))
allocate(start_lon(nbergs_in_file))
allocate(start_lat(nbergs_in_file))
allocate(start_day(nbergs_in_file))
allocate(start_mass(nbergs_in_file))
allocate(mass_scaling(nbergs_in_file))
allocate(mass_of_bits(nbergs_in_file))
allocate(static_berg(nbergs_in_file))
allocate(heat_density(nbergs_in_file))
allocate(ine(nbergs_in_file))
allocate(jne(nbergs_in_file))
allocate(start_year(nbergs_in_file))
if (replace_iceberg_num) then
call error_mesg('read_restart_bergs', "Calculating new iceberg ID's", WARNING)
allocate(iceberg_num(nbergs_in_file))
else
allocate(id_cnt(nbergs_in_file))
allocate(id_ij(nbergs_in_file))
endif
call read_unlimited_axis(filename,'lon',lon,domain=grd%domain)
call read_unlimited_axis(filename,'lat',lat,domain=grd%domain)
call read_unlimited_axis(filename,'uvel',uvel,domain=grd%domain)
call read_unlimited_axis(filename,'vvel',vvel,domain=grd%domain)
call read_unlimited_axis(filename,'mass',mass,domain=grd%domain)
call read_real_vector(filename,'axn',axn,grd%domain,value_if_not_in_file=0.)
call read_real_vector(filename,'ayn',ayn,grd%domain,value_if_not_in_file=0.)
call read_real_vector(filename,'bxn',bxn,grd%domain,value_if_not_in_file=0.)
call read_real_vector(filename,'byn',byn,grd%domain,value_if_not_in_file=0.)
call read_unlimited_axis(filename,'thickness',thickness,domain=grd%domain)
call read_unlimited_axis(filename,'width',width,domain=grd%domain)
call read_unlimited_axis(filename,'length',length,domain=grd%domain)
call read_unlimited_axis(filename,'start_lon',start_lon,domain=grd%domain)
call read_unlimited_axis(filename,'start_lat',start_lat,domain=grd%domain)
call read_unlimited_axis(filename,'start_day',start_day,domain=grd%domain)
call read_unlimited_axis(filename,'start_mass',start_mass,domain=grd%domain)
call read_unlimited_axis(filename,'mass_scaling',mass_scaling,domain=grd%domain)
call read_real_vector(filename,'mass_of_bits',mass_of_bits,domain=grd%domain,value_if_not_in_file=0.)
call read_real_vector(filename,'heat_density',heat_density,domain=grd%domain,value_if_not_in_file=0.)
call read_unlimited_axis(filename,'ine',ine,domain=grd%domain)
call read_unlimited_axis(filename,'jne',jne,domain=grd%domain)
call read_unlimited_axis(filename,'start_year',start_year,domain=grd%domain)
if (replace_iceberg_num) then
call read_int_vector(filename,'iceberg_num',iceberg_num,grd%domain,value_if_not_in_file=-1)
else
call read_int_vector(filename,'id_cnt',id_cnt,grd%domain)
call read_int_vector(filename,'id_ij',id_ij,grd%domain)
endif
call read_real_vector(filename,'static_berg',static_berg,grd%domain,value_if_not_in_file=0.)
endif
! Find approx outer bounds for tile
lon0=minval( grd%lon(grd%isc-1:grd%iec,grd%jsc-1:grd%jec) )
lon1=maxval( grd%lon(grd%isc-1:grd%iec,grd%jsc-1:grd%jec) )
lat0=minval( grd%lat(grd%isc-1:grd%iec,grd%jsc-1:grd%jec) )
lat1=maxval( grd%lat(grd%isc-1:grd%iec,grd%jsc-1:grd%jec) )
do k=1, nbergs_in_file
localberg%lon=lon(k)
localberg%lat=lat(k)
if (.not. ignore_ij_restart) then ! read i,j position and avoid the "find" step
localberg%ine=ine(k)
localberg%jne=jne(k)
if ( localberg%ine>=grd%isc .and. localberg%ine<=grd%iec .and. &
localberg%jne>=grd%jsc .and.localberg%jne<=grd%jec ) then
lres=.true.
else
lres=.false.
endif
else ! i,j are not available from the file so we search the grid to find out if we reside on this PE
if (use_slow_find) then
lres=find_cell(grd, localberg%lon, localberg%lat, localberg%ine, localberg%jne)
else
lres=find_cell_by_search(grd, localberg%lon, localberg%lat, localberg%ine, localberg%jne)
endif
endif
if (really_debug) then
write(stderrunit,'(a,i8,a,2f9.4,a,i8)') 'diamonds, read_restart_bergs: berg ',k,' is at ',localberg%lon,localberg%lat,&
& ' on PE ',mpp_pe()
write(stderrunit,*) 'diamonds, read_restart_bergs: lres = ',lres
endif
if (lres) then ! true if we reside on this PE grid
localberg%uvel=uvel(k)
localberg%vvel=vvel(k)
localberg%mass=mass(k)
localberg%axn=axn(k) !Alon
localberg%ayn=ayn(k) !Alon
localberg%uvel_old=uvel(k) !Alon
localberg%vvel_old=vvel(k) !Alon
localberg%lon_old=lon(k) !Alon
localberg%lat_old=lat(k) !Alon
localberg%bxn=bxn(k) !Alon
localberg%byn=byn(k) !Alon
localberg%thickness=thickness(k)
localberg%width=width(k)
localberg%length=length(k)
localberg%start_lon=start_lon(k)
localberg%start_lat=start_lat(k)
localberg%start_year=start_year(k)
if (replace_iceberg_num) then
localberg%id = generate_id(grd, localberg%ine, localberg%jne)
else
localberg%id=id_from_2_ints(id_cnt(k), id_ij(k))
endif
localberg%start_day=start_day(k)
localberg%start_mass=start_mass(k)
localberg%mass_scaling=mass_scaling(k)
localberg%mass_of_bits=mass_of_bits(k)
localberg%halo_berg=0.
localberg%static_berg=static_berg(k)
localberg%heat_density=heat_density(k)
localberg%first_bond=>null()
if (really_debug) lres=is_point_in_cell(grd, localberg%lon, localberg%lat, localberg%ine, localberg%jne, explain=.true.)
lres=pos_within_cell(grd, localberg%lon, localberg%lat, localberg%ine, localberg%jne, localberg%xi, localberg%yj)
!call add_new_berg_to_list(bergs%first, localberg, quick=.true.)
if (bergs%grd%area(localberg%ine,localberg%jne) .ne. 0) then
call add_new_berg_to_list(bergs%list(localberg%ine,localberg%jne)%first, localberg)
else
if (mpp_pe().eq.mpp_root_pe()) then
print * , 'Grounded iceberg: ', lat(k),lon(k), localberg%id
call error_mesg('diamonds, read_restart_bergs', 'Iceberg not added because it is grounded', WARNING)
endif
endif
if (really_debug) call print_berg(stderrunit, bergs%list(localberg%ine,localberg%jne)%first, 'read_restart_bergs, add_new_berg_to_list')
elseif (multiPErestart .and. io_tile_id(1) .lt. 0) then
call error_mesg('diamonds, read_restart_bergs', 'berg in PE file was not on PE!', FATAL)
endif
enddo
if (found_restart) then
deallocate(lon, &
lat, &
uvel, &
vvel, &
mass, &
axn, &
ayn, &
bxn, &
byn, &
thickness, &
width, &
length, &
start_lon, &
start_lat, &
start_day, &
start_mass, &
mass_scaling, &
mass_of_bits, &
static_berg, &
heat_density, &
ine, &
jne, &
start_year )
if (replace_iceberg_num) then
deallocate(iceberg_num)
else
deallocate(id_cnt)
deallocate(id_ij)
endif
! This block only works for IO_LAYOUT=1,1 or 0,0 but not for arbitrary layouts.
! I'm commenting this out until we find a way to implement the same sorts of checks
! that work for all i/o layouts. -AJA
!Checking the total number of icebergs read from the restart file.
!nbergs_read=count_bergs(bergs)
!call mpp_sum(nbergs_read)
!if (mpp_pe().eq.mpp_root_pe()) then
! write(*,'(a,i8,a,i8,a)') 'diamonds, read_restart_bergs: Number of Icebergs in restart file=',nbergs_in_file,' Number of Icebergs read=', nbergs_read
! if (nbergs_read .gt. nbergs_in_file) then
! call error_mesg('diamonds, read_restart_bergs', 'More icebergs read than exist in restart file.', FATAL)
! elseif (nbergs_read .lt. nbergs_in_file) then
! if (bergs%ignore_missing_restart_bergs) then
! call error_mesg('diamonds, read_restart_bergs', 'Some Icebergs from restart file were not found (ignore_missing flag is on)', WARNING)
! else
! call error_mesg('diamonds, read_restart_bergs', 'Some Icebergs from restart file were not found', FATAL)
! endif
! elseif (nbergs_read .eq. nbergs_in_file) then
! write(*,'(a,i8,a,i8,a)') 'diamonds, read_restart_bergs: Number of icebergs read (#',nbergs_read,') matches the number of icebergs in the file'
! endif
!endif
elseif(.not. found_restart .and. bergs%nbergs_start==0 .and. generate_test_icebergs) then
call generate_bergs(bergs,Time)
endif
call check_for_duplicates_in_parallel(bergs)
bergs%floating_mass_start=sum_mass(bergs)
call mpp_sum( bergs%floating_mass_start )
bergs%icebergs_mass_start=sum_mass(bergs,justbergs=.true.)
call mpp_sum( bergs%icebergs_mass_start )
bergs%bergy_mass_start=sum_mass(bergs,justbits=.true.)
call mpp_sum( bergs%bergy_mass_start )
if (mpp_pe().eq.mpp_root_pe().and.verbose) write(*,'(a)') 'diamonds, read_restart_bergs: completed'
end subroutine read_restart_bergs
!> Read a vector of reals from file and use a default value if variable is missing
subroutine read_real_vector(filename, varname, values, domain, value_if_not_in_file)
character(len=*), intent(in) :: filename !< Name of file to read from
character(len=*), intent(in) :: varname !< Name of variable to read
real, intent(out) :: values(:) !< Returned vector of reals
type(domain2D), intent(in) :: domain !< Parallel decomposition
real, optional, intent(in) :: value_if_not_in_file !< Value to use if variable is not in file
if (present(value_if_not_in_file).and..not.field_exist(filename, varname)) then
values(:)=value_if_not_in_file
else
call read_unlimited_axis(filename,varname,values,domain=domain)
endif
end subroutine read_real_vector
!> Read a vector of integers from file and use a default value if variable is missing
subroutine read_int_vector(filename, varname, values, domain, value_if_not_in_file)
character(len=*), intent(in) :: filename !< Name of file to read from
character(len=*), intent(in) :: varname !< Name of variable to read
integer, intent(out) :: values(:) !< Returned vector of integers
type(domain2D), intent(in) :: domain !< Parallel decomposition
integer, optional, intent(in) :: value_if_not_in_file !< Value to use if variable is not in file
if (present(value_if_not_in_file).and..not.field_exist(filename, varname)) then
values(:)=value_if_not_in_file
else
call read_unlimited_axis(filename,varname,values,domain=domain)
endif
end subroutine read_int_vector
!> Generate bergs for the purpose of debugging
subroutine generate_bergs(bergs,Time)
! Arguments
type(icebergs), pointer :: bergs !< Icebergs container
type(time_type), intent(in) :: Time !< Model time
! Local variables
type(icebergs_gridded), pointer :: grd
integer :: i,j
type(iceberg) :: localberg ! NOT a pointer but an actual local variable
integer :: iyr, imon, iday, ihr, imin, isec
logical :: lres
! For convenience
grd=>bergs%grd
call get_date(Time, iyr, imon, iday, ihr, imin, isec)
do j=grd%jsc,grd%jec; do i=grd%isc,grd%iec
if (grd%msk(i,j)>0. .and. abs(grd%latc(i,j))>80.0) then
if (max(grd%lat(i,j),grd%lat(i-1,j),grd%lat(i,j-1),grd%lat(i-1,j-1))>89.999) cycle ! Cannot use this at Pole cells
localberg%xi=-999.
localberg%yj=-999.
localberg%ine=i
localberg%jne=j
localberg%lon=grd%lonc(i,j)
localberg%lat=grd%latc(i,j)
localberg%lon_old=localberg%lon
localberg%lat_old=localberg%lat
localberg%mass=bergs%initial_mass(1)
localberg%thickness=bergs%initial_thickness(1)
localberg%width=bergs%initial_width(1)
localberg%length=bergs%initial_length(1)
localberg%start_lon=localberg%lon
localberg%start_lat=localberg%lat
localberg%start_year=iyr
localberg%start_day=float(iday)+(float(ihr)+float(imin)/60.)/24.
localberg%start_mass=localberg%mass
localberg%mass_scaling=bergs%mass_scaling(1)
localberg%mass_of_bits=0.
localberg%halo_berg=0.
localberg%static_berg=0.
localberg%heat_density=0.
localberg%axn=0. !Alon
localberg%ayn=0. !Alon
localberg%uvel_old=0. !Alon
localberg%vvel_old=0. !Alon
localberg%bxn=0. !Alon
localberg%byn=0. !Alon
!Berg A
call loc_set_berg_pos(grd, 0.9, 0.5, 1., 0., localberg)
localberg%id = generate_id(grd, i, j)
call add_new_berg_to_list(bergs%list(i,j)%first, localberg)
!Berg B
call loc_set_berg_pos(grd, 0.1, 0.5, -1., 0., localberg)
localberg%id = generate_id(grd, i, j)
call add_new_berg_to_list(bergs%list(i,j)%first, localberg)
!Berg C
call loc_set_berg_pos(grd, 0.5, 0.9, 0., 1., localberg)
localberg%id = generate_id(grd, i, j)
call add_new_berg_to_list(bergs%list(i,j)%first, localberg)
!Berg D
call loc_set_berg_pos(grd, 0.5, 0.1, 0., -1., localberg)
localberg%id = generate_id(grd, i, j)
call add_new_berg_to_list(bergs%list(i,j)%first, localberg)
endif
enddo; enddo
bergs%nbergs_start=count_bergs(bergs)
call mpp_sum(bergs%nbergs_start)
if (mpp_pe().eq.mpp_root_pe()) &
write(*,'(a,i8,a)') 'diamonds, generate_bergs: ',bergs%nbergs_start,' were generated'
end subroutine generate_bergs
subroutine loc_set_berg_pos(grd, xi, yj, uvel, vvel, berg)
type(icebergs_gridded), pointer :: grd !< Container for gridded fields
real, intent(in) :: xi !< Non-dimensional x-position within cell to give berg
real, intent(in) :: yj !< Non-dimensional y-position within cell to give berg
real, intent(in) :: uvel !< Zonal velocity to give berg
real, intent(in) :: vvel !< Meridional velocity to give berg
type(iceberg), intent(inout) :: berg !< An iceberg
integer :: i, j
logical :: lres
i = berg%ine ; j = berg%jne
if (max(grd%lat(i,j),grd%lat(i-1,j),grd%lat(i,j-1),grd%lat(i-1,j-1))>89.999) then
berg%lon=grd%lonc(i,j)
berg%lat=grd%latc(i,j)
berg%xi=0.5 ; berg%yj=0.5
else
berg%lon=bilin(grd, grd%lon, i, j, xi, yj)
berg%lat=bilin(grd, grd%lat, i, j, xi, yj)
berg%xi=xi ; berg%yj=yj
endif
berg%uvel=uvel ; berg%vvel=vvel
berg%lon_old=berg%lon ; berg%lat_old=berg%lat
berg%start_lon=berg%lon ; berg%start_lat=berg%lat
lres=pos_within_cell(grd, berg%lon, berg%lat, berg%ine, berg%jne, berg%xi, berg%yj)
if (.not. lres) then
lres=pos_within_cell(grd, berg%lon, berg%lat, berg%ine, berg%jne, berg%xi, berg%yj, explain=.true.)
write(0,*) lres, i, j, xi, yj, uvel, vvel
write(0,*) lres, berg%ine, berg%jne, berg%xi, berg%yj
write(0,*) 'bx=',berg%lon, 'gx=',grd%lon(i-1,j-1), grd%lon(i,j-1), grd%lon(i,j), grd%lon(i-1,j),'cx=', grd%lonc(i,j)
write(0,*) 'by=',berg%lat, 'gy=',grd%lat(i-1,j-1), grd%lat(i,j-1), grd%lat(i,j), grd%lat(i-1,j),'cy=', grd%latc(i,j)
stop 'generate_bergs, loc_set_berg_pos(): VERY FATAL!'
endif
end subroutine loc_set_berg_pos
!> Read bond restart file
subroutine read_restart_bonds(bergs,Time)
! Arguments
type(icebergs), pointer :: bergs !< Icebergs container
type(time_type), intent(in) :: Time !< Model time
! Local variables
integer :: k, siz(4), nbonds_in_file
logical :: lres, found_restart, found
logical :: first_berg_found, second_berg_found
logical :: multiPErestart ! Not needed with new restart read; currently kept for compatibility
real :: lon0, lon1, lat0, lat1
character(len=33) :: filename, filename_base
type(icebergs_gridded), pointer :: grd
type(iceberg) :: localberg ! NOT a pointer but an actual local variable
type(iceberg) , pointer :: this, first_berg, second_berg
type(bond) , pointer :: current_bond
integer :: stderrunit
integer :: number_first_bonds_matched !How many first bond bergs found on pe
integer :: number_second_bonds_matched !How many second bond bergs found on pe
integer :: number_perfect_bonds ! How many complete bonds formed
integer :: number_partial_bonds ! How many either complete/partial bonds formed.
integer :: number_perfect_bonds_with_first_on_pe ! How many bonds with first bond on the compuational domain
integer :: all_pe_number_perfect_bonds, all_pe_number_partial_bonds
integer :: all_pe_number_first_bonds_matched, all_pe_number_second_bonds_matched
integer :: all_pe_number_perfect_bonds_with_first_on_pe
integer :: ine, jne
logical :: search_data_domain
real :: berg_found, berg_found_all_pe
integer, allocatable, dimension(:) :: id_cnt, id_ij, &
first_berg_jne, &
first_berg_ine, &
other_berg_jne, &
other_berg_ine
integer(kind=8), allocatable, dimension(:) :: first_id, &
other_id
!integer, allocatable, dimension(:,:) :: iceberg_counter_grd
! Get the stderr unit number
stderrunit=stderr()
! For convenience
grd=>bergs%grd
! Zero out nbergs_in_file
nbonds_in_file = 0
all_pe_number_perfect_bonds=0
filename_base=trim(restart_input_dir)//'bonds_iceberg.res.nc'
found_restart = find_restart_file(filename_base, filename, multiPErestart, io_tile_id(1))
call error_mesg('read_restart_bonds_bergs_new', 'Using icebergs bond restart read', NOTE)
filename = filename_base
call get_field_size(filename,'i',siz, field_found=found, domain=bergs%grd%domain)
nbonds_in_file = siz(1)
if (mpp_pe() .eq. mpp_root_pe()) then
write(stderrunit,*) 'diamonds, bond read restart : ','Number of bonds in file', nbonds_in_file
endif
if (nbonds_in_file .gt. 0) then
allocate(first_id(nbonds_in_file))
allocate(other_id(nbonds_in_file))
allocate(id_cnt(nbonds_in_file))
allocate(id_ij(nbonds_in_file))
allocate(first_berg_jne(nbonds_in_file))
allocate(first_berg_ine(nbonds_in_file))
allocate(other_berg_ine(nbonds_in_file))
allocate(other_berg_jne(nbonds_in_file))
call read_unlimited_axis(filename,'first_id_cnt',id_cnt,domain=grd%domain)
call read_unlimited_axis(filename,'first_id_ij',id_ij,domain=grd%domain)
do k=1, nbonds_in_file
first_id(k) = id_from_2_ints( id_cnt(k), id_ij(k) )
enddo
call read_unlimited_axis(filename,'other_id_cnt',id_cnt,domain=grd%domain)
call read_unlimited_axis(filename,'other_id_ij',id_ij,domain=grd%domain)
do k=1, nbonds_in_file
other_id(k) = id_from_2_ints( id_cnt(k), id_ij(k) )
enddo
deallocate(id_cnt, id_ij)
call read_unlimited_axis(filename,'first_berg_jne',first_berg_jne,domain=grd%domain)
call read_unlimited_axis(filename,'first_berg_ine',first_berg_ine,domain=grd%domain)
call read_unlimited_axis(filename,'other_berg_jne',other_berg_jne,domain=grd%domain)
call read_unlimited_axis(filename,'other_berg_ine',other_berg_ine,domain=grd%domain)
number_first_bonds_matched=0
number_second_bonds_matched=0
number_perfect_bonds=0
number_partial_bonds=0
number_perfect_bonds_with_first_on_pe=0
do k=1, nbonds_in_file
! If i,j in restart files are not good, then we find the berg position of the bond addresses manually:
if (ignore_ij_restart) then
!Finding first iceberg in bond
ine=999 ; jne=999 ; berg_found=0.0 ; search_data_domain=.true.
call find_individual_iceberg(bergs,first_id(k), ine, jne,berg_found,search_data_domain)
berg_found_all_pe=berg_found
call mpp_sum(berg_found_all_pe)
if (berg_found_all_pe .gt. 0.5) then
first_berg_ine(k)=ine
first_berg_jne(k)=jne
else
print * , 'First bond berg not located: ', first_id(k),berg_found, mpp_pe(),ine, jne
call error_mesg('read_restart_bonds_bergs_new', 'First iceberg in bond not found on any pe', FATAL)
endif
!else
!Finding other iceberg other iceberg
ine=999 ; jne=999 ; berg_found=0.0 ; search_data_domain =.true.
call find_individual_iceberg(bergs,other_id(k), ine, jne, berg_found,search_data_domain)
berg_found_all_pe=berg_found
call mpp_sum(berg_found_all_pe)
if (berg_found_all_pe .gt. 0.5) then
!if (berg_found_all_pe .gt. 1.5) then
! call error_mesg('read_restart_bonds_bergs_new', 'Other iceberg bond found on more than one pe', FATAL)
!else
other_berg_ine(k)=ine
other_berg_jne(k)=jne
!endif
else
call error_mesg('read_restart_bonds_bergs_new', 'Other iceberg in bond not found on any pe', FATAL)
endif
if (berg_found_all_pe .lt. 0.5) then
print * , 'First bond berg not located: ', other_id(k),berg_found, mpp_pe(),ine, jne
call error_mesg('read_restart_bonds_bergs_new', 'First bond iceberg not located', FATAL)
endif
endif
! Decide whether the first iceberg is on the processeor