forked from NOAA-GFDL/SIS2
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathSIS_dyn_cgrid.F90
1751 lines (1581 loc) · 88.9 KB
/
SIS_dyn_cgrid.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
module SIS_dyn_cgrid
!***********************************************************************
!* 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 *
!***********************************************************************
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
! !
! C-grid SEA ICE DYNAMICS using ELASTIC-VISCOUS-PLASTIC RHEOLOGY adapted from !
! Hunke and Dukowicz (JPO 1997, H&D hereafter) with some derivation from SIS1 !
! and with guidance from the C-grid implementations of sea-ice in MITgcm as !
! documented in MITgcm user notes by Martin Losch and in LIM3 by S. Bouillon !
! et al. (Ocean Modelling, 2009 & 2013). This code initially written by !
! Robert Hallberg in 2013. !
! !
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
use SIS_diag_mediator, only : post_SIS_data, SIS_diag_ctrl
use SIS_diag_mediator, only : query_SIS_averaging_enabled, enable_SIS_averaging
use SIS_diag_mediator, only : register_diag_field=>register_SIS_diag_field
use SIS_debugging, only : chksum, Bchksum, hchksum, vec_chksum_C
use SIS_debugging, only : check_redundant_B, check_redundant_C, vec_chksum_C
use MOM_error_handler, only : SIS_error=>MOM_error, FATAL, WARNING, NOTE, SIS_mesg=>MOM_mesg
use MOM_file_parser, only : get_param, log_param, read_param, log_version, param_file_type
use MOM_domains, only : pass_var, pass_vector, CGRID_NE, CORNER, pe_here
use MOM_hor_index, only : hor_index_type
use MOM_io, only : open_file
use MOM_io, only : APPEND_FILE, ASCII_FILE, MULTIPLE, SINGLE_FILE
use SIS_hor_grid, only : SIS_hor_grid_type
use fms_io_mod, only : register_restart_field, restart_file_type
use mpp_domains_mod, only : domain2D
use time_manager_mod, only : time_type, set_time, operator(+), operator(-)
use time_manager_mod, only : set_date, get_time, get_date
implicit none ; private
#include <SIS2_memory.h>
public :: SIS_C_dyn_init, SIS_C_dynamics, SIS_C_dyn_end, SIS_C_dyn_register_restarts
type, public :: SIS_C_dyn_CS ; private
real, allocatable, dimension(:,:) :: &
str_t, & ! The tension stress tensor component, in Pa m.
str_d, & ! The divergence stress tensor component, in Pa m.
str_s ! The shearing stress tensor component (cross term), in Pa m.
! parameters for calculating water drag and internal ice stresses
logical :: SLAB_ICE = .false. ! should we do old style GFDL slab ice?
real :: p0 = 2.75e4 ! pressure constant (Pa)
real :: p0_rho ! The pressure constant divided by ice density, N m kg-1.
real :: c0 = 20.0 ! another pressure constant
real :: cdw = 3.24e-3 ! ice/water drag coef. (nondim)
real :: EC = 2.0 ! yield curve axis ratio
real :: Rho_ocean = 1030.0 ! The nominal density of sea water, in kg m-3.
real :: Rho_ice = 905.0 ! The nominal density of sea ice, in kg m-3.
real :: drag_bg_vel2 = 0.0 ! A background (subgridscale) velocity for drag
! with the ocean squared, in m2 s-2.
real :: min_ocn_inertial_h = 0. ! A minimum ocean thickness used to limit the viscous coupling
! rate implied for the ocean by the ice-ocean stress.
real :: Tdamp ! The damping timescale of the stress tensor components
! toward their equilibrium solution due to the elastic terms,
! in s.
real :: del_sh_min_scale = 2.0 ! A scaling factor for the minimum permitted
! value of minimum shears used in the denominator
! of the stress equations, nondim. I suspect that
! this needs to be greater than 1.
real :: CFL_trunc ! Velocity components will be truncated when they
! are large enough that the corresponding CFL number
! exceeds this value, nondim.
logical :: CFL_check_its ! If true, check the CFL number for every iteration
! of the rheology solver; otherwise only check the
! final velocities that are used for transport.
logical :: specified_ice ! If true, the sea ice is specified and there is
! no need for ice dynamics.
logical :: debug ! If true, write verbose checksums for debugging purposes.
logical :: debug_redundant ! If true, debug redundant points.
logical :: project_drag_vel ! If true, project forward the ice velocity used
! in the drag calculation to avoid an instability
! that can occur when an finite stress is applied
! to thin ice moving with the velocity of the ocean.
logical :: project_ci ! If true, project the ice concentration and
! related ice strength changes due to the convergent
! or divergent ice flow.
logical :: weak_coast_stress = .false.
logical :: weak_low_shear = .false.
integer :: evp_sub_steps ! The number of iterations in the EVP dynamics
! for each slow time step.
real :: dt_Rheo ! The maximum sub-cycling time step for the EVP dynamics.
type(time_type), pointer :: Time ! A pointer to the ice model's clock.
type(SIS_diag_ctrl), pointer :: diag ! A structure that is used to regulate the
! timing of diagnostic output.
integer, pointer :: ntrunc ! The number of times the velocity has been truncated
! since the last call to write_ice_statistics.
character(len = 200) :: u_trunc_file ! The complete path to files in which a
character(len = 200) :: v_trunc_file ! column's worth of accelerations are
! written if velocity truncations occur.
integer :: u_file, v_file ! The unit numbers for opened u- or v- truncation
! files, or -1 if they have not yet been opened.
integer :: cols_written ! The number of columns whose output has been
! written by this PE during the current run.
integer :: max_writes ! The maximum number of times any PE can write out
! a column's worth of accelerations during a run.
logical :: FirstCall = .true.
integer :: id_fix = -1, id_fiy = -1, id_fcx = -1, id_fcy = -1
integer :: id_fwx = -1, id_fwy = -1, id_sigi = -1, id_sigii = -1
integer :: id_stren = -1, id_stren0 = -1
integer :: id_ui = -1, id_vi = -1, id_Coru = -1, id_Corv = -1
integer :: id_PFu = -1, id_PFv = -1, id_fpx = -1, id_fpy = -1
integer :: id_fix_d = -1, id_fix_t = -1, id_fix_s = -1
integer :: id_fiy_d = -1, id_fiy_t = -1, id_fiy_s = -1
integer :: id_str_d = -1, id_str_t = -1, id_str_s = -1
integer :: id_sh_d = -1, id_sh_t = -1, id_sh_s = -1
integer :: id_del_sh = -1, id_del_sh_min = -1
integer :: id_mis = -1, id_ci = -1, id_ci0 = -1, id_miu = -1, id_miv = -1
integer :: id_ui_hifreq = -1, id_vi_hifreq = -1
integer :: id_str_d_hifreq = -1, id_str_t_hifreq = -1, id_str_s_hifreq = -1
integer :: id_sh_d_hifreq = -1, id_sh_t_hifreq = -1, id_sh_s_hifreq = -1
integer :: id_sigi_hifreq = -1, id_sigii_hifreq = -1
integer :: id_stren_hifreq = -1, id_ci_hifreq = -1
integer :: id_siu = -1, id_siv = -1, id_sispeed = -1 ! SIMIP diagnostics
end type SIS_C_dyn_CS
contains
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
! SIS_C_dyn_init - initialize the ice dynamics and set parameters. !
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
subroutine SIS_C_dyn_init(Time, G, param_file, diag, CS, ntrunc)
type(time_type), target, intent(in) :: Time
type(SIS_hor_grid_type), intent(in) :: G
type(param_file_type), intent(in) :: param_file
type(SIS_diag_ctrl), target, intent(inout) :: diag
type(SIS_C_dyn_CS), pointer :: CS
integer, target, optional, intent(inout) :: ntrunc
! Arguments: Time - The current model time.
! (in) G - The ocean's grid structure.
! (in) param_file - A structure indicating the open file to parse for
! model parameter values.
! (in) diag - A structure that is used to regulate diagnostic output.
! (in/out) CS - A pointer that is set to point to the control structure
! for this module.
! (in/out,opt) ntrunc - The integer that stores the number of times the velocity
! has been truncated since the last call to write_ice_statistics.
! This subroutine sets the parameters and registers the diagnostics associated
! with the ice dynamics.
! This include declares and sets the variable "version".
#include "version_variable.h"
character(len=40) :: mod = "SIS_C_dyn" ! This module's name.
real, parameter :: missing = -1e34
if (.not.associated(CS)) then
call SIS_error(FATAL, "SIS_C_dyn_init called with an unassociated control structure. \n"//&
"SIS_C_dyn_register_restarts must be called before SIS_C_dyn_init.")
return
endif
CS%diag => diag
CS%Time => Time
if (present(ntrunc)) then ; CS%ntrunc => ntrunc ; else ; allocate(CS%ntrunc) ; endif
CS%ntrunc = 0
! Read all relevant parameters and write them to the model log.
call log_version(param_file, mod, version)
call get_param(param_file, mod, "SPECIFIED_ICE", CS%specified_ice, &
"If true, the ice is specified and there is no dynamics.", &
default=.false.)
if ( CS%specified_ice ) then
CS%evp_sub_steps = 0 ; CS%dt_Rheo = -1.0
call log_param(param_file, mod, "NSTEPS_DYN", CS%evp_sub_steps, &
"The number of iterations in the EVP dynamics for each \n"//&
"slow time step. With SPECIFIED_ICE this is always 0.")
else
call get_param(param_file, mod, "DT_RHEOLOGY", CS%dt_Rheo, &
"The sub-cycling time step for iterating the rheology \n"//&
"and ice momentum equations. If DT_RHEOLOGY is negative, \n"//&
"the time step is set via NSTEPS_DYN.", units="seconds", &
default=-1.0)
CS%evp_sub_steps = -1
if (CS%dt_Rheo <= 0.0) &
call get_param(param_file, mod, "NSTEPS_DYN", CS%evp_sub_steps, &
"The number of iterations of the rheology and ice \n"//&
"momentum equations for each slow ice time step.", default=432)
endif
call get_param(param_file, mod, "ICE_TDAMP_ELASTIC", CS%Tdamp, &
"The damping timescale associated with the elastic terms \n"//&
"in the sea-ice dynamics equations (if positive) or the \n"//&
"fraction of DT_ICE_DYNAMICS (if negative).", &
units = "s or nondim", default=-0.2)
call get_param(param_file, mod, "WEAK_LOW_SHEAR_ICE", CS%weak_low_shear, &
"If true, the divergent stresses go toward 0 in the C-grid \n"//&
"dynamics when the shear magnitudes are very weak. \n"//&
"Otherwise they go to -P_ice. This setting is temporary.", &
default=.false.)
call get_param(param_file, mod, "PROJECT_ICE_DRAG_VEL", CS%project_drag_vel, &
"If true, project forward the ice velocity used in the \n"//&
"drag calculation to avoid an instability that can occur \n"//&
"when an finite stress is applied to thin ice moving with \n"//&
"the velocity of the ocean.", default=.true.)
call get_param(param_file, mod, "ICE_YIELD_ELLIPTICITY", CS%EC, &
"The ellipticity coefficient for the plastic yield curve \n"//&
"in the sea-ice rheology. For an infinite ellipticity \n"//&
"(i.e., a cavitating fluid rheology), use 0.", &
units = "Nondim", default=2.0)
call get_param(param_file, mod, "ICE_STRENGTH_PSTAR", CS%p0, &
"A constant in the expression for the ice strength, \n"//&
"P* in Hunke & Dukowicz 1997.", units="Pa", default=2.75e4)
call get_param(param_file, mod, "ICE_STRENGTH_CSTAR", CS%c0, &
"A constant in the exponent of the expression for the \n"//&
"ice strength, c* in Hunke & Dukowicz 1997.", &
units="nondim", default=20.)
call get_param(param_file, mod, "ICE_CDRAG_WATER", CS%cdw, &
"The drag coefficient between the sea ice and water.", &
units="nondim", default=3.24e-3)
call get_param(param_file, mod, "MIN_OCN_INTERTIAL_H", CS%min_ocn_inertial_h, &
"A minimum ocean thickness used to limit the viscous coupling rate\n"//&
"implied for the ocean by the ice-ocean stress. Only used if positive.", &
units="m", default=0.0)
call get_param(param_file, mod, "ICE_DEL_SH_MIN_SCALE", CS%del_sh_min_scale, &
"A scaling factor for the lower bound on the shear rates \n"//&
"used in the denominator of the stress calculation. This \n"//&
"probably needs to be greater than 1.", units="nondim", default=2.0)
call get_param(param_file, mod, "PROJECT_ICE_CONCENTRATION", CS%project_ci, &
"If true, project the evolution of the ice concentration \n"//&
"due to the convergence or divergence of the ice flow.", default=.true.)
call get_param(param_file, mod, "RHO_OCEAN", CS%Rho_ocean, &
"The nominal density of sea water as used by SIS.", &
units="kg m-3", default=1030.0)
call get_param(param_file, mod, "RHO_ICE", CS%Rho_ice, &
"The nominal density of sea ice as used by SIS.", &
units="kg m-3", default=905.0)
CS%p0_rho = CS%p0 / CS%Rho_ice
call get_param(param_file, mod, "CFL_TRUNCATE", CS%CFL_trunc, &
"The value of the CFL number that will cause ice velocity \n"//&
"components to be truncated; instability can occur past 0.5.", &
units="nondim", default=0.5)
call get_param(param_file, mod, "CFL_TRUNC_DYN_ITS", CS%CFL_check_its, &
"If true, check the CFL number for every iteration of the \n"//&
"rheology solver; otherwise only the final velocities that \n"//&
"are used for transport are checked.", &
default=.false.)
call get_param(param_file, mod, "DEBUG", CS%debug, &
"If true, write out verbose debugging data.", default=.false.)
call get_param(param_file, mod, "DEBUG_REDUNDANT", CS%debug_redundant, &
"If true, debug redundant data points.", default=CS%debug)
if ( CS%specified_ice ) then
CS%slab_ice = .true.
call log_param(param_file, mod, "USE_SLAB_ICE", CS%slab_ice, &
"Use the very old slab-style ice. With SPECIFIED_ICE, \n"//&
"USE_SLAB_ICE is always true.")
else
call get_param(param_file, mod, "USE_SLAB_ICE", CS%slab_ice, &
"If true, use the very old slab-style ice.", default=.false.)
endif
call get_param(param_file, mod, "U_TRUNC_FILE", CS%u_trunc_file, &
"The absolute path to the file where the accelerations \n"//&
"leading to zonal velocity truncations are written. \n"//&
"Leave this empty for efficiency if this diagnostic is \n"//&
"not needed.", default="")
call get_param(param_file, mod, "V_TRUNC_FILE", CS%v_trunc_file, &
"The absolute path to the file where the accelerations \n"//&
"leading to meridional velocity truncations are written. \n"//&
"Leave this empty for efficiency if this diagnostic is \n"//&
"not needed.", default="")
call get_param(param_file, mod, "MAX_TRUNC_FILE_SIZE_PER_PE", CS%max_writes, &
"The maximum number of colums of truncations that any PE \n"//&
"will write out during a run.", default=50)
! if (len_trim(dirs%output_directory) > 0) then
! if (len_trim(CS%u_trunc_file) > 0) &
! CS%u_trunc_file = trim(dirs%output_directory)//trim(CS%u_trunc_file)
! if (len_trim(CS%v_trunc_file) > 0) &
! CS%v_trunc_file = trim(dirs%output_directory)//trim(CS%v_trunc_file)
! call log_param(param_file, mod, "output_dir/U_TRUNC_FILE", CS%u_trunc_file)
! call log_param(param_file, mod, "output_dir/V_TRUNC_FILE", CS%v_trunc_file)
! endif
CS%u_file = -1 ; CS%v_file = -1 ; CS%cols_written = 0
CS%id_sigi = register_diag_field('ice_model','SIGI' ,diag%axesT1, Time, &
'first stress invariant', 'none', missing_value=missing)
CS%id_sigii = register_diag_field('ice_model','SIGII' ,diag%axesT1, Time, &
'second stress invariant', 'none', missing_value=missing)
CS%id_stren = register_diag_field('ice_model','STRENGTH' ,diag%axesT1, Time, &
'ice strength', 'Pa*m', missing_value=missing)
CS%id_stren0 = register_diag_field('ice_model','STREN_0' ,diag%axesT1, Time, &
'ice strength at start of rheology', 'Pa*m', missing_value=missing)
CS%id_fix = register_diag_field('ice_model', 'FI_X', diag%axesCu1, Time, &
'ice internal stress - x component', 'Pa', missing_value=missing, &
interp_method='none')
CS%id_fiy = register_diag_field('ice_model', 'FI_Y', diag%axesCv1, Time, &
'ice internal stress - y component', 'Pa', missing_value=missing, &
interp_method='none')
CS%id_fcx = register_diag_field('ice_model', 'FC_X', diag%axesCu1, Time, &
'Coriolis force - x component', 'Pa', missing_value=missing, &
interp_method='none')
CS%id_fcy = register_diag_field('ice_model', 'FC_Y', diag%axesCv1, Time, &
'Coriolis force - y component', 'Pa', missing_value=missing, &
interp_method='none')
CS%id_Coru = register_diag_field('ice_model', 'Cor_ui', diag%axesCu1, Time,&
'Coriolis ice acceleration - x component', 'm s-2', &
missing_value=missing, interp_method='none')
CS%id_Corv = register_diag_field('ice_model', 'Cor_vi', diag%axesCv1, Time,&
'Coriolis ice acceleration - y component', 'm s-2', &
missing_value=missing, interp_method='none')
CS%id_fpx = register_diag_field('ice_model', 'FP_X', diag%axesCu1, Time, &
'Pressure force - x component', 'Pa', missing_value=missing, &
interp_method='none')
CS%id_fpy = register_diag_field('ice_model', 'FP_Y', diag%axesCv1, Time, &
'Pressure force - y component', 'Pa', missing_value=missing, &
interp_method='none')
CS%id_PFu = register_diag_field('ice_model', 'Pfa_ui', diag%axesCu1, Time, &
'Pressure-force ice acceleration - x component', 'm s-2', &
missing_value=missing, interp_method='none')
CS%id_PFv = register_diag_field('ice_model', 'Pfa_vi', diag%axesCv1, Time, &
'Pressure-force ice acceleration - y component', 'm s-2', &
missing_value=missing, interp_method='none')
CS%id_fwx = register_diag_field('ice_model', 'FW_X', diag%axesCu1, Time, &
'water stress on ice - x component', 'Pa', missing_value=missing, &
interp_method='none')
CS%id_fwy = register_diag_field('ice_model', 'FW_Y', diag%axesCv1, Time, &
'water stress on ice - y component', 'Pa', missing_value=missing, &
interp_method='none')
CS%id_ui = register_diag_field('ice_model', 'UI', diag%axesCu1, Time, &
'ice velocity - x component', 'm/s', missing_value=missing, &
interp_method='none')
CS%id_vi = register_diag_field('ice_model', 'VI', diag%axesCv1, Time, &
'ice velocity - y component', 'm/s', missing_value=missing, &
interp_method='none')
CS%id_mis = register_diag_field('ice_model', 'MIS_tot', diag%axesT1, Time, &
'Mass of ice and snow at t-points', 'kg m-2', missing_value=missing)
CS%id_ci0 = register_diag_field('ice_model', 'CI_tot', diag%axesT1, Time, &
'Initial summed concentration of ice at t-points', 'nondim', &
missing_value=missing)
CS%id_ci = register_diag_field('ice_model', 'CI_proj', diag%axesT1, Time, &
'Projected summed concentration of ice at t-points', 'nondim', &
missing_value=missing)
CS%id_miu = register_diag_field('ice_model', 'MI_U', diag%axesCu1, Time, &
'Mass of ice and snow at u-points', 'kg m-2', &
missing_value=missing, interp_method='none')
CS%id_miv = register_diag_field('ice_model', 'MI_V', diag%axesCv1, Time, &
'Mass of ice and snow at v-points', 'kg m-2', &
missing_value=missing, interp_method='none')
CS%id_fix_d = register_diag_field('ice_model', 'FI_d_X', diag%axesCu1, Time, &
'ice divergence internal stress - x component', 'Pa', missing_value=missing, &
interp_method='none')
CS%id_fiy_d = register_diag_field('ice_model', 'FI_d_Y', diag%axesCv1, Time, &
'ice divergence internal stress - y component', 'Pa', missing_value=missing, &
interp_method='none')
CS%id_fix_t = register_diag_field('ice_model', 'FI_t_X', diag%axesCu1, Time, &
'ice tension internal stress - x component', 'Pa', missing_value=missing, &
interp_method='none')
CS%id_fiy_t = register_diag_field('ice_model', 'FI_t_Y', diag%axesCv1, Time, &
'ice tension internal stress - y component', 'Pa', missing_value=missing, &
interp_method='none')
CS%id_fix_s = register_diag_field('ice_model', 'FI_s_X', diag%axesCu1, Time, &
'ice shearing internal stress - x component', 'Pa', missing_value=missing, &
interp_method='none')
CS%id_fiy_s = register_diag_field('ice_model', 'FI_s_Y', diag%axesCv1, Time, &
'ice shearing internal stress - y component', 'Pa', missing_value=missing, &
interp_method='none')
CS%id_str_d = register_diag_field('ice_model', 'str_d', diag%axesT1, Time, &
'ice divergence internal stress', 'Pa', missing_value=missing)
CS%id_str_t = register_diag_field('ice_model', 'str_t', diag%axesT1, Time, &
'ice tension internal stress', 'Pa', missing_value=missing)
CS%id_str_s = register_diag_field('ice_model', 'str_s', diag%axesB1, Time, &
'ice shearing internal stress', 'Pa', missing_value=missing)
CS%id_sh_d = register_diag_field('ice_model', 'sh_d', diag%axesT1, Time, &
'ice divergence strain rate', 's-1', missing_value=missing)
CS%id_sh_t = register_diag_field('ice_model', 'sh_t', diag%axesT1, Time, &
'ice tension strain rate', 's-1', missing_value=missing)
CS%id_sh_s = register_diag_field('ice_model', 'sh_s', diag%axesB1, Time, &
'ice shearing strain rate', 's-1', missing_value=missing)
CS%id_del_sh = register_diag_field('ice_model', 'del_sh', diag%axesT1, Time, &
'ice strain rate magnitude', 's-1', missing_value=missing)
CS%id_del_sh_min = register_diag_field('ice_model', 'del_sh_min', diag%axesT1, Time, &
'minimum ice strain rate magnitude', 's-1', missing_value=missing)
CS%id_ui_hifreq = register_diag_field('ice_model', 'ui_hf', diag%axesCu1, Time, &
'ice velocity - x component', 'm/s', missing_value=missing, &
interp_method='none')
CS%id_vi_hifreq = register_diag_field('ice_model', 'vi_hf', diag%axesCv1, Time, &
'ice velocity - y component', 'm/s', missing_value=missing, &
interp_method='none')
CS%id_str_d_hifreq = register_diag_field('ice_model', 'str_d_hf', diag%axesT1, Time, &
'ice divergence internal stress', 'Pa', missing_value=missing)
CS%id_str_t_hifreq = register_diag_field('ice_model', 'str_t_hf', diag%axesT1, Time, &
'ice tension internal stress', 'Pa', missing_value=missing)
CS%id_str_s_hifreq = register_diag_field('ice_model', 'str_s_hf', diag%axesB1, Time, &
'ice shearing internal stress', 'Pa', missing_value=missing)
CS%id_sh_d_hifreq = register_diag_field('ice_model', 'sh_d_hf', diag%axesT1, Time, &
'ice divergence rate', 's-1', missing_value=missing)
CS%id_sh_t_hifreq = register_diag_field('ice_model', 'sh_t_hf', diag%axesT1, Time, &
'ice tension rate', 's-1', missing_value=missing)
CS%id_sh_s_hifreq = register_diag_field('ice_model', 'sh_s_hf', diag%axesB1, Time, &
'ice shearing rate', 's-1', missing_value=missing)
CS%id_sigi_hifreq = register_diag_field('ice_model','sigI_hf' ,diag%axesT1, Time, &
'first stress invariant', 'none', missing_value=missing)
CS%id_sigii_hifreq = register_diag_field('ice_model','sigII_hf' ,diag%axesT1, Time, &
'second stress invariant', 'none', missing_value=missing)
CS%id_ci_hifreq = register_diag_field('ice_model', 'CI_hf', diag%axesT1, Time, &
'Summed concentration of ice at t-points', 'nondim', missing_value=missing)
CS%id_stren_hifreq = register_diag_field('ice_model','STRENGTH_hf' ,diag%axesT1, Time, &
'ice strength', 'Pa*m', missing_value=missing)
CS%id_siu = register_diag_field('ice_model', 'siu', diag%axesT1, Time, &
'ice velocity - x component', 'm/s', missing_value=missing, &
interp_method='none')
CS%id_siv = register_diag_field('ice_model', 'siv', diag%axesT1, Time, &
'ice velocity - y component', 'm/s', missing_value=missing, &
interp_method='none')
CS%id_sispeed = register_diag_field('ice_model', 'sispeed', diag%axesT1, Time, &
'ice speed', 'm/s', missing_value=missing)
end subroutine SIS_C_dyn_init
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
! find_ice_strength - magnitude of force on ice in plastic deformation !
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
subroutine find_ice_strength(mi, ci, ice_strength, G, CS, halo_sz) ! ??? may change to do loop
type(SIS_hor_grid_type), intent(in) :: G
real, dimension(SZI_(G),SZJ_(G)), intent(in) :: mi, ci
real, dimension(SZI_(G),SZJ_(G)), intent(out) :: ice_strength
type(SIS_C_dyn_CS), pointer :: CS
integer, optional, intent(in) :: halo_sz
integer :: i, j, isc, iec, jsc, jec, halo
halo = 0 ; if (present(halo_sz)) halo = halo_sz
isc = G%isc-halo ; iec = G%iec+halo ; jsc = G%jsc-halo ; jec = G%jec+halo
do j=jsc,jec ; do i=isc,iec
ice_strength(i,j) = CS%p0_rho*mi(i,j)*exp(-CS%c0*max(1.0-ci(i,j),0.0))
enddo ; enddo
end subroutine find_ice_strength
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
! SIS_C_dynamics - take a single dynamics timestep with EVP subcycles !
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
subroutine SIS_C_dynamics(ci, msnow, mice, ui, vi, uo, vo, &
fxat, fyat, sea_lev, fxoc, fyoc, dt_slow, G, CS)
type(SIS_hor_grid_type), intent(inout) :: G
real, dimension(SZI_(G),SZJ_(G)), intent(in ) :: ci, msnow, mice ! ice properties
real, dimension(SZIB_(G),SZJ_(G)), intent(inout) :: ui ! ice velocity
real, dimension(SZI_(G),SZJB_(G)), intent(inout) :: vi ! ice velocity
real, dimension(SZIB_(G),SZJ_(G)), intent(in ) :: uo ! ocean velocity
real, dimension(SZI_(G),SZJB_(G)), intent(in ) :: vo ! ocean velocity
real, dimension(SZIB_(G),SZJ_(G)), intent(in ) :: fxat ! air stress on ice
real, dimension(SZI_(G),SZJB_(G)), intent(in ) :: fyat ! air stress on ice
real, dimension(SZI_(G),SZJ_(G)), intent(in ) :: sea_lev ! sea level
real, dimension(SZIB_(G),SZJ_(G)), intent( out) :: fxoc ! ice stress on ocean
real, dimension(SZI_(G),SZJB_(G)), intent( out) :: fyoc ! ice stress on ocean
real, intent(in ) :: dt_slow
type(SIS_C_dyn_CS), pointer :: CS
! Arguments: ci - The sea ice concentration, nondim.
! (in) msnow - The mass per unit total area (ice covered and ice free)
! of the snow, in kg m-2.
! (in) mice - The mass per unit total area (ice covered and ice free)
! of the ice, in kg m-2.
! (inout) ui - The zonal ice velocity, in m s-1.
! (inout) vi - The meridional ice velocity, in m s-1.
! (in) uo - The zonal ocean velocity, in m s-1.
! (in) vo - The meridional ocean velocity, in m s-1.
! (in) sea_lev - The height of the sea level, including contributions
! from non-levitating ice from an earlier time step, in m.
! (in) dt_slow - The amount of time over which the ice dynamics are to be
! advanced, in s.
! (in) G - The ocean's grid structure.
! (in/out) CS - A pointer to the control structure for this module.
real, dimension(SZI_(G),SZJ_(G)) :: &
sh_Dt, & ! sh_Dt is the horizontal tension (du/dx - dv/dy) including
! all metric terms, in s-1.
sh_Dd ! sh_Dd is the flow divergence (du/dx + dv/dy) including all
! metric terms, in s-1.
real, dimension(SZIB_(G),SZJB_(G)) :: &
sh_Ds ! sh_Ds is the horizontal shearing strain (du/dy + dv/dx)
! including all metric terms, in s-1.
real, dimension(SZI_(G),SZJ_(G)) :: &
mis, & ! Total snow and ice mass per unit area, in kg m-2.
pres_mice, & ! The ice internal pressure per unit column mass, in N m / kg.
ci_proj, & ! The projected ice concentration, nondim.
zeta, & ! The ice bulk viscosity, in Pa m s or N s / m.
del_sh, & ! The magnitude of the shear rates, in s-1.
diag_val, & ! A temporary diagnostic array.
del_sh_min_pr, & ! When multiplied by pres_mice, this gives the minimum
! value of del_sh that is used in the calculation of zeta,
! in s-1. This is set based on considerations of numerical
! stability, and varies with the grid spacing.
dx2T, dy2T, & ! dx^2 or dy^2 at T points, in m2.
dx_dyT, dy_dxT, & ! dx/dy or dy_dx at T points, nondim.
siu, siv, sispeed ! diagnostics on T points, m/s
real, dimension(SZIB_(G),SZJ_(G)) :: &
fxic, & ! Zonal force due to internal stresses, in Pa.
fxic_d, fxic_t, fxic_s, &
ui_min_trunc, & ! The range of v-velocities beyond which the velocities
ui_max_trunc, & ! are truncated, in m s-1, or 0 for land cells.
Cor_u, & ! Zonal Coriolis acceleration, in m s-2.
PFu, & ! Zonal hydrostatic pressure driven acceleration, in m s-2.
u_tmp, & ! A temporary copy of the old values of ui, in m s-1.
u_IC, & ! The initial zonal ice velocities, in m s-1.
mi_u, & ! The total ice and snow mass interpolated to u points, in kg m-2.
f2dt_u, &! The squared effective Coriolis parameter at u-points times a
! time step, in s-1.
I1_f2dt2_u ! 1 / ( 1 + f^2 dt^2) at u-points, nondimensional.
real, dimension(SZI_(G),SZJB_(G)) :: &
fyic, & ! Meridional force due to internal stresses, in Pa.
fyic_d, fyic_t, fyic_s, &
vi_min_trunc, & ! The range of v-velocities beyond which the velocities
vi_max_trunc, & ! are truncated, in m s-1, or 0 for land cells.
Cor_v, & ! Meridional Coriolis acceleration, in m s-2.
PFv, & ! Meridional hydrostatic pressure driven acceleration, in m s-2.
v_IC, & ! The initial meridional ice velocities, in m s-1.
mi_v, & ! The total ice and snow mass interpolated to v points, in kg m-2.
f2dt_v, &! The squared effective Coriolis parameter at v-points times a
! time step, in s-1.
I1_f2dt2_v ! 1 / ( 1 + f^2 dt^2) at v-points, nondimensional.
real, dimension(SZIB_(G),SZJB_(G)) :: &
mi_ratio_A_q, & ! A ratio of the masses interpolated to the faces around a
! vorticity point that ranges between (4 mi_min/mi_max) and 1,
! divided by the sum of the ocean areas around a point, in m-2.
q, & ! A potential-vorticity-like field for the ice, the Coriolis
! parameter divided by a spatially averaged mass per unit area,
! in s-1 m2 kg-1.
dx2B, dy2B, & ! dx^2 or dy^2 at B points, in m2.
dx_dyB, dy_dxB ! dx/dy or dy_dx at B points, nondim.
real, dimension(SZIB_(G),SZJ_(G)) :: &
azon, bzon, & ! _zon & _mer are the values of the Coriolis force which
czon, dzon, & ! are applied to the neighboring values of vi & ui,
amer, bmer, & ! respectively to get the barotropic inertial rotation,
cmer, dmer ! in units of s-1. azon and amer couple the same pair of
! velocities, but with the influence going in opposite
! directions.
real :: Cor ! A Coriolis accleration, in m s-2.
real :: fxic_now, fyic_now ! ice internal stress convergence, in kg m-1 s-2.
real :: drag_u, drag_v ! Drag rates with the ocean at u & v points, in kg m-2 s-1.
real :: drag_max ! A maximum drag rate allowed in the ocean, in kg m-2 s-1.
real :: tot_area ! The sum of the area of the four neighboring cells, in m2.
real :: dxharm ! The harmonic mean of the x- and y- grid spacings, in m.
real :: muq2, mvq2 ! The product of the u- and v-face masses per unit cell
! area surrounding a vorticity point, in kg2 m-4.
real :: muq, mvq ! The u- and v-face masses per unit cell area extrapolated
! to a vorticity point on the coast, in kg m-2.
real :: pres_sum ! The sum of the internal ice pressures aroung a point, in Pa.
real :: min_rescale ! The smallest of the 4 surrounding values of rescale, ND.
real :: I_1pdt_T ! 1.0 / (1.0 + dt_2Tdamp)
real :: I_1pE2dt_T ! 1.0 / (1.0 + EC^2 * dt_2Tdamp)
real :: v2_at_u ! The squared v-velocity interpolated to u points, in m s-1.
real :: u2_at_v ! The squared u-velocity interpolated to v points, in m s-1.
real :: uio_init, m_uio_explicit, uio_pred ! , uio
real :: vio_init, m_vio_explicit, vio_pred ! , vio
real :: I_cdRhoDt, cdRho
real :: b_vel0 ! The initial difference between the velocity magnitude
! and the absolute value of the u- or v- component, plus
! the ice thickness divided by the time step and the drag
! coefficient, all in m s-1.
real :: uio_C ! A u-velocity difference between the ocean and ice, in m s-1.
real :: vio_C ! A v-velocity difference between the ocean and ice, in m s-1.
real :: Tdamp ! The damping timescale of the stress tensor components
! toward their equilibrium solution due to the elastic terms,
! in s.
real :: dt ! The short timestep associated with the EVP dynamics, in s.
real :: dt_2Tdamp ! The ratio of the timestep to the elastic damping timescale.
real :: dt_cumulative ! The elapsed time within this call to EVP dynamics, in s.
integer :: EVP_steps ! The number of EVP sub-steps that will actually be taken.
real :: I_sub_steps ! The number inverse of the number of EVP time steps per
! slow time step.
real :: EC2 ! EC^2, where EC is the yield curve axis ratio.
real :: I_EC2 ! 1/EC^2, where EC is the yield curve axis ratio.
real :: I_EC ! 1/EC, where EC is the yield curve axis ratio.
real :: I_2EC ! 1/(2*EC), where EC is the yield curve axis ratio.
real, parameter :: H_subroundoff = 1e-30 ! A negligible thickness, in m, that
! can be cubed without underflow.
real :: m_neglect ! A tiny mass per unit area, in kg m-2.
real :: m_neglect2 ! A tiny mass per unit area squared, in kg2 m-4.
real :: m_neglect4 ! A tiny mass per unit area to the 4th power, in kg4 m-8.
real :: sum_area ! The sum of ocean areas around a vorticity point, in m2.
type(time_type) :: &
time_it_start, & ! The starting time of the iteratve steps.
time_step_end, & ! The end time of an iterative step.
time_end_in ! The end time for diagnostics when this routine started.
real :: time_int_in ! The diagnostics' time interval when this routine started.
logical :: do_hifreq_output ! If true, output occurs every iterative step.
logical :: do_trunc_its ! If true, overly large velocities in the iterations are truncated.
integer :: halo_sh_Ds ! The halo size that can be used in calculating sh_Ds.
integer :: i, j, isc, iec, jsc, jec, n
isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec
if (.not.associated(CS)) call SIS_error(FATAL, &
"SIS_C_dynamics: Module must be initialized before it is used.")
if ((isc - G%isdB < 2) .or. (jsc - G%jsdB < 2)) call SIS_error(FATAL, &
"SIS_C_dynamics is written to require a 2-point halo or 1-point and symmetric memory.")
halo_sh_Ds = min(isc-G%isd, jsc-G%jsd, 2)
! Zero these arrays to accumulate sums.
fxoc(:,:) = 0.0 ; fyoc(:,:) = 0.0
fxic(:,:) = 0.0 ; fyic(:,:) = 0.0
Cor_u(:,:) = 0.0 ; Cor_v(:,:) = 0.0
fxic_d(:,:) = 0.0 ; fyic_d(:,:) = 0.0 ; fxic_t(:,:) = 0.0 ; fyic_t(:,:) = 0.0
fxic_s(:,:) = 0.0 ; fyic_s(:,:) = 0.0
if (CS%SLAB_ICE) then
ui(:,:) = uo(:,:) ; vi(:,:) = vo(:,:)
fxoc(:,:) = fxat(:,:) ; fyoc(:,:) = fyat(:,:)
return
end if
if ((CS%evp_sub_steps<=0) .and. (CS%dt_Rheo<=0.0)) return
if (CS%FirstCall) then
! These halo updates are only needed if the str_... arrays have just
! been read from a restart file. Otherwise the halos contain valid data.
call pass_var(CS%str_d, G%Domain) ; call pass_var(CS%str_t, G%Domain)
call pass_var(CS%str_s, G%Domain, position=CORNER)
CS%FirstCall = .false.
endif
if (CS%dt_Rheo > 0.0) then
EVP_steps = max(CEILING(dt_slow/CS%dt_Rheo - 0.0001), 1)
else
EVP_steps = CS%evp_sub_steps
endif
dt = dt_slow/EVP_steps
drag_max = CS%Rho_ocean * CS%min_ocn_inertial_h / dt_slow
I_cdRhoDt = 1.0 / (CS%cdw * CS%Rho_ocean * dt)
do_trunc_its = (CS%CFL_check_its .and. (CS%CFL_trunc > 0.0) .and. (dt_slow > 0.0))
EC2 = CS%EC**2
I_EC = 0.0 ; if (CS%EC > 0.0) I_EC = 1.0 / CS%EC
I_2EC = 0.0 ; if (CS%EC > 0.0) I_2EC = 0.5 / CS%EC
I_EC2 = 0.0 ; if (EC2 > 0.0) I_EC2 = 1.0 / EC2
do_hifreq_output = .false.
if ((CS%id_ui_hifreq > 0) .or. (CS%id_vi_hifreq > 0) .or. &
(CS%id_str_d_hifreq > 0) .or. (CS%id_str_t_hifreq > 0) .or. &
(CS%id_str_s_hifreq > 0) .or. (CS%id_sh_d_hifreq > 0) .or. &
(CS%id_sh_t_hifreq > 0) .or. (CS%id_sh_s_hifreq > 0) .or. &
(CS%id_ci_hifreq > 0) .or. (CS%id_stren_hifreq > 0)) then
do_hifreq_output = query_SIS_averaging_enabled(CS%diag, time_int_in, time_end_in)
if (do_hifreq_output) &
time_it_start = time_end_in - set_time(int(floor(dt_slow+0.5)))
endif
Tdamp = CS%Tdamp
if (CS%Tdamp == 0.0) then
! Hunke (2001) chooses a specified multiple (0.36) of dt_slow for Tdamp,
! and shows that stability requires Tdamp > 2*dt. Here 0.2 is used instead
! for greater stability.
Tdamp = max(0.2*dt_slow, 3.0*dt)
elseif (CS%Tdamp < 0.0) then
Tdamp = max(-CS%Tdamp*dt_slow, 3.0*dt)
endif
dt_2Tdamp = dt / (2.0 * Tdamp)
ui_min_trunc(:,:) = 0.0 ; ui_max_trunc(:,:) = 0.0
vi_min_trunc(:,:) = 0.0 ; vi_max_trunc(:,:) = 0.0
m_neglect = H_subroundoff*CS%Rho_ice
m_neglect2 = m_neglect**2 ; m_neglect4 = m_neglect**4
!$OMP parallel default(none) shared(isc,iec,jsc,jec,G,CS,dt_slow,ui_min_trunc,u_IC,ui, &
!$OMP ui_max_trunc,vi_min_trunc,vi_max_trunc,v_IC,vi,mice, &
!$OMP msnow,ci,dt,Tdamp,I_2EC,mis,ci_proj,pres_mice, &
!$OMP dx2B,dy2B,dx_dyB,dy_dxB,dx2T,dy2T,dx_dyT,dy_dxT, &
!$OMP mi_ratio_A_q,m_neglect4,m_neglect2,mi_u,mi_v,q, &
!$OMP m_neglect,azon,bzon,czon,dzon,f2dt_u,I1_f2dt2_u,PFu, &
!$OMP sea_lev,amer,bmer,cmer,dmer,f2dt_v,I1_f2dt2_v,PFv, &
!$OMP del_sh_min_pr ) &
!$OMP private(dxharm,sum_area,muq2,mvq2,muq,mvq,tot_area)
if ((CS%CFL_trunc > 0.0) .and. (dt_slow > 0.0)) then
!$OMP do
do j=jsc,jec
do I=isc-1,iec ; if (G%dy_Cu(I,j) > 0.0) then
ui_min_trunc(I,j) = (-CS%CFL_trunc) * G%areaT(i+1,j) / (dt_slow*G%dy_Cu(I,j))
ui_max_trunc(I,j) = CS%CFL_trunc * G%areaT(i,j) / (dt_slow*G%dy_Cu(I,j))
endif ; enddo
do I=isc-1,iec ; u_IC(I,j) = ui(I,j) ; enddo
enddo
!$OMP end do nowait
!$OMP do
do J=jsc-1,jec
do i=isc,iec ; if (G%dx_Cv(i,J) > 0.0) then
vi_min_trunc(i,J) = (-CS%CFL_trunc) * G%areaT(i,j+1) / (dt_slow*G%dx_Cv(i,J))
vi_max_trunc(i,J) = CS%CFL_trunc * G%areaT(i,j) / (dt_slow*G%dx_Cv(i,J))
endif ; enddo
do i=isc,iec ; v_IC(i,J) = vi(i,j) ; enddo
enddo
!$OMP end do nowait
endif
!$OMP do
do j=jsc-1,jec+1 ; do i=isc-1,iec+1
! Store the total snow and ice mass.
mis(i,j) = mice(i,j) + msnow(i,j)
ci_proj(i,j) = ci(i,j)
! Precompute pres_mice and the minimum value of del_sh for stability.
pres_mice(i,j) = CS%p0_rho*exp(-CS%c0*max(1.0-ci(i,j),0.0))
dxharm = 2.0*G%dxT(i,j)*G%dyT(i,j) / (G%dxT(i,j) + G%dyT(i,j))
! Setting a minimum value of del_sh is sufficient to guarantee numerical
! stability of the overall time-stepping for the velocities and stresses.
! Setting a minimum value of the shear magnitudes is equivalent to setting
! a maximum value of the effective lateral viscosities.
! I think that this is stable when CS%del_sh_min_scale >= 1. -RWH
if (dxharm > 0.) then
del_sh_min_pr(i,j) = (2.0*CS%del_sh_min_scale * dt**2) / (Tdamp * dxharm**2)
else
del_sh_min_pr(i,j) = 0.
endif
enddo ; enddo
! Ensure that the input stresses are not larger than could be justified by
! the ice pressure now, as the ice might have melted or been advected away
! during the thermodynamic and transport phases.
call limit_stresses(pres_mice, mice, CS%str_d, CS%str_t, CS%str_s, G, CS)
! Zero out ice velocities with no mass.
!$OMP do
do j=jsc,jec ; do I=isc-1,iec
if (G%mask2dCu(I,j) * (mis(i,j)+mis(i+1,j)) == 0.0) ui(I,j) = 0.0
enddo ; enddo
!$OMP end do nowait
!$OMP do
do J=jsc-1,jec ; do i=isc,iec
if (G%mask2dCv(i,J) * (mis(i,j)+mis(i,j+1)) == 0.0) vi(I,j) = 0.0
enddo ; enddo
!$OMP end do nowait
!$OMP do
do J=jsc-1,jec ; do I=isc-1,iec
dx2B(I,J) = G%dxBu(I,J)*G%dxBu(I,J) ; dy2B(I,J) = G%dyBu(I,J)*G%dyBu(I,J)
enddo ; enddo
!$OMP end do nowait
!$OMP do
do J=jsc-2,jec+1 ; do I=isc-2,iec+1
dx_dyB(I,J) = G%dxBu(I,J)*G%IdyBu(I,J) ; dy_dxB(I,J) = G%dyBu(I,J)*G%IdxBu(I,J)
enddo ; enddo
!$OMP end do nowait
!$OMP do
do j=jsc-1,jec+1 ; do i=isc-1,iec+1
dx2T(i,j) = G%dxT(i,j)*G%dxT(i,j) ; dy2T(i,j) = G%dyT(i,j)*G%dyT(i,j)
dx_dyT(i,j) = G%dxT(i,j)*G%IdyT(i,j) ; dy_dxT(i,j) = G%dyT(i,j)*G%IdxT(i,j)
enddo ; enddo
!$OMP end do nowait
!$OMP do
do J=jsc-1,jec ; do I=isc-1,iec
if (CS%weak_coast_stress) then
sum_area = (G%areaT(i,j) + G%areaT(i+1,j+1)) + (G%areaT(i,j+1) + G%areaT(i+1,j))
else
sum_area = (G%mask2dT(i,j)*G%areaT(i,j) + G%mask2dT(i+1,j+1)*G%areaT(i+1,j+1)) + &
(G%mask2dT(i,j+1)*G%areaT(i,j+1) + G%mask2dT(i+1,j)*G%areaT(i+1,j))
endif
if (sum_area <= 0.0) then
! This is a land point.
mi_ratio_A_q(I,J) = 0.0
elseif (G%mask2dBu(I,J)>0.5) then
! This is an interior ocean point.
! Determine an appropriately averaged mass on q-points. The following
! expression for mi_q is mi when the masses are all equal, and goes to 4
! times the smallest mass averaged onto the 4 adjacent velocity points. It
! comes from taking the harmonic means of the harmonic means of the
! arithmetic mean masses at the velocity points. mi_ratio goes from 4 times
! the ratio of the smallest mass at velocity points over the largest mass
! at velocity points up to 1.
muq2 = 0.25 * (mis(i,j) + mis(i+1,j)) * (mis(i,j+1) + mis(i+1,j+1))
mvq2 = 0.25 * (mis(i,j) + mis(i,j+1)) * (mis(i+1,j) + mis(i+1,j+1))
mi_ratio_A_q(I,J) = 32.0 * muq2 * mvq2 / ((m_neglect4 + (muq2 + mvq2) * &
((mis(i,j) + mis(i+1,j+1)) + (mis(i,j+1) + mis(i+1,j)))**2) * &
sum_area)
elseif ((G%mask2dCu(I,j) + G%mask2dCu(I,j+1)) + &
(G%mask2dCv(i,J) + G%mask2dCv(i+1,J)) > 1.5) then
! This is a corner point, and there are 1 unmasked u-point and 1 v-point.
! The ratio below goes from 4 times the ratio of the smaller of the two
! masses at velocity points over the larger up to 1.
muq = 0.5 * (G%mask2dCu(I,j) * (mis(i,j) + mis(i+1,j)) + &
G%mask2dCu(I,j+1) * (mis(i,j+1) + mis(i+1,j+1)) )
mvq = 0.5 * (G%mask2dCv(i,J) * (mis(i,j) + mis(i,j+1)) + &
G%mask2dCv(i+1,J) * (mis(i+1,j) + mis(i+1,j+1)) )
mi_ratio_A_q(I,J) = 4.0 * muq * mvq / ((m_neglect2 + (muq + mvq)**2) * sum_area)
else
! This is a straight coastline or all neighboring velocity points are
! masked out. In any case, with just 1 point, the ratio is always 1.
mi_ratio_A_q(I,J) = 1.0 / sum_area
endif
enddo ; enddo
!$OMP end do nowait
!$OMP do
do j=jsc-1,jec+1 ; do I=isc-1,iec
mi_u(I,j) = 0.5*(mis(i+1,j) + mis(i,j))
enddo ; enddo
!$OMP end do nowait
!$OMP do
do J=jsc-1,jec ; do i=isc-1,iec+1
mi_v(i,J) = 0.5*(mis(i,j+1) + mis(i,j))
enddo ; enddo
!$OMP end do nowait
!$OMP do
do J=jsc-1,jec ; do I=isc-1,iec
tot_area = (G%areaT(i,j) + G%areaT(i+1,j+1)) + (G%areaT(i+1,j) + G%areaT(i,j+1))
q(I,J) = G%CoriolisBu(I,J) * tot_area / &
(((G%areaT(i,j) * mis(i,j) + G%areaT(i+1,j+1) * mis(i+1,j+1)) + &
(G%areaT(i+1,j) * mis(i+1,j) + G%areaT(i,j+1) * mis(i,j+1))) + tot_area * m_neglect)
enddo ; enddo
!$OMP do
do j=jsc,jec ; do I=isc-1,iec
! Calculate terms related to the Coriolis force on the zonal velocity.
azon(I,j) = 0.25 * mi_v(i+1,J) * q(I,J)
bzon(I,j) = 0.25 * mi_v(i,J) * q(I,J)
czon(I,j) = 0.25 * mi_v(i,J-1) * q(I,J-1)
dzon(I,j) = 0.25 * mi_v(i+1,J-1) * q(I,J-1)
f2dt_u(I,j) = dt * 4.0 * ((azon(I,j)**2 + czon(I,j)**2) + &
(bzon(I,j)**2 + dzon(I,j)**2))
I1_f2dt2_u(I,j) = 1.0 / ( 1.0 + dt * f2dt_u(I,j) )
! Calculate the zonal acceleration due to the sea level slope.
PFu(I,j) = -G%g_Earth*(sea_lev(i+1,j)-sea_lev(i,j)) * G%IdxCu(I,j)
enddo ; enddo
!$OMP end do nowait
!$OMP do
do J=jsc-1,jec ; do i=isc,iec
! Calculate terms related to the Coriolis force on the meridional velocity.
amer(I-1,j) = 0.25 * mi_u(I-1,j) * q(I-1,J)
bmer(I,j) = 0.25 * mi_u(I,j) * q(I,J)
cmer(I,j+1) = 0.25 * mi_u(I,j+1) * q(I,J)
dmer(I-1,j+1) = 0.25 * mi_u(I-1,j+1) * q(I-1,J)
f2dt_v(i,J) = dt * 4.0 * ((amer(I-1,j)**2 + cmer(I,j+1)**2) + &
(bmer(I,j)**2 + dmer(I-1,j+1)**2))
I1_f2dt2_v(i,J) = 1.0 / ( 1.0 + dt * f2dt_v(i,J) )
! Calculate the meridional acceleration due to the sea level slope.
PFv(i,J) = -G%g_Earth*(sea_lev(i,j+1)-sea_lev(i,j)) * G%IdyCv(i,J)
enddo ; enddo
!$OMP end parallel
if (CS%debug .or. CS%debug_redundant) then
call vec_chksum_C("PF[uv] in SIS_C_dynamics", PFu, PFv, G)
call vec_chksum_C("f[xy]at in SIS_C_dynamics", fxat, fyat, G)
call vec_chksum_C("[uv]i pre-steps SIS_C_dynamics", ui, vi, G)
call vec_chksum_C("[uv]o in SIS_C_dynamics", uo, vo, G)
endif
dt_cumulative = 0.0
! Do the iterative time steps.
do n=1,EVP_steps
dt_cumulative = dt_cumulative + dt
! If there is a 2-point wide halo and symmetric memory, this is the only
! halo update that is needed per iteration. With a 1-point wide halo and
! symmetric memory, an update is also needed for sh_Ds.
call pass_vector(ui, vi, G%Domain, stagger=CGRID_NE)
! Calculate the strain tensor for viscosities and forcing elastic eqn.
! The following are the forms of the horizontal tension and hori-
! shearing strain advocated by Smagorinsky (1993) and discussed in
! Griffies and Hallberg (MWR, 2000). Similar forms are used in the sea
! ice model of Bouillon et al. (Ocean Modelling, 2009).
! The calculation of sh_Ds has the widest halo. The logic below avoids
! a halo update when possible.
! With a halo of >= 2 this is: do J=jsc-2,jec+1 ; do I=isc-2,iec+1
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,halo_sh_Ds,sh_Ds,G, &
!$OMP dx_dyB,dy_dxB,ui,vi)
do J=jsc-halo_sh_Ds,jec+halo_sh_Ds-1 ; do I=isc-halo_sh_Ds,iec+halo_sh_Ds-1
! This uses a no-slip boundary condition.
sh_Ds(I,J) = (2.0-G%mask2dBu(I,J)) * &
(dx_dyB(I,J)*(ui(I,j+1)*G%IdxCu(I,j+1) - ui(I,j)*G%IdxCu(I,j)) + &
dy_dxB(I,J)*(vi(i+1,J)*G%IdyCv(i+1,J) - vi(i,J)*G%IdyCv(i,J)))
enddo ; enddo
if (halo_sh_Ds < 2) call pass_var(sh_Ds, G%Domain, position=CORNER)
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,sh_Dt,sh_Dd,dy_dxT,dx_dyT,G,ui,vi)
do j=jsc-1,jec+1 ; do i=isc-1,iec+1
sh_Dt(i,j) = (dy_dxT(i,j)*(G%IdyCu(I,j) * ui(I,j) - &
G%IdyCu(I-1,j)*ui(I-1,j)) - &
dx_dyT(i,j)*(G%IdxCv(i,J) * vi(i,J) - &
G%IdxCv(i,J-1)*vi(i,J-1)))
sh_Dd(i,j) = (G%IareaT(i,j)*(G%dyCu(I,j) * ui(I,j) - &
G%dyCu(I-1,j)*ui(I-1,j)) + &
G%IareaT(i,j)*(G%dxCv(i,J) * vi(i,J) - &
G%dxCv(i,J-1)*vi(i,J-1)))
enddo ; enddo
if (CS%project_ci) then
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,ci_proj,ci,dt_cumulative, &
!$OMP sh_Dd,pres_mice,CS)
do j=jsc-1,jec+1 ; do i=isc-1,iec+1
! Estimate future ice concentrations from the approximate expression
! d ci / dt = - ci * sh_Dt
! The choice to base this on the final velocity, the initial concentration
! and the elapsed time is because it is that final velocity that will drive
! ice convergence.
ci_proj(i,j) = ci(i,j) * exp(-dt_cumulative*sh_Dd(i,j))
! Recompute pres_mice.
pres_mice(i,j) = CS%p0_rho*exp(-CS%c0*max(1.0-ci_proj(i,j),0.0))
enddo ; enddo
endif
! calculate viscosities - how often should we do this ?
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,del_sh,zeta,sh_Dd,sh_Dt, &
!$OMP I_EC2,sh_Ds,pres_mice,mice,del_sh_min_pr)
do j=jsc-1,jec+1 ; do i=isc-1,iec+1
! Averaging the squared shearing strain is larger than squaring
! the averaged strain. I don't know what is better. -RWH
del_sh(i,j) = sqrt(sh_Dd(i,j)**2 + I_EC2 * (sh_Dt(i,j)**2 + &
(0.25 * ((sh_Ds(I-1,J-1) + sh_Ds(I,J)) + &
(sh_Ds(I-1,J) + sh_Ds(I,J-1))))**2 ) ) ! H&D eqn 9
if (max(del_sh(i,j), del_sh_min_pr(i,j)*pres_mice(i,j)) /=0.) then
zeta(i,j) = 0.5*pres_mice(i,j)*mice(i,j) / &
max(del_sh(i,j), del_sh_min_pr(i,j)*pres_mice(i,j))
else
zeta(i,j) = 0.
endif
enddo ; enddo
! Step the stress component equations semi-implicitly.
I_1pdt_T = 1.0 / (1.0 + dt_2Tdamp)
I_1pE2dt_T = 1.0 / (1.0 + EC2*dt_2Tdamp)
if (CS%weak_low_shear) then
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,CS,I_1pdt_T,dt_2Tdamp,zeta, &
!$OMP sh_Dd,del_sh,I_EC2,sh_Dt)
do j=jsc-1,jec+1 ; do i=isc-1,iec+1
! This expression uses that Pres=2*del_sh*zeta with an elliptic yield curve.
CS%str_d(i,j) = I_1pdt_T * ( CS%str_d(i,j) + dt_2Tdamp * &
( zeta(i,j) * (sh_Dd(i,j) - del_sh(i,j)) ) )
CS%str_t(i,j) = I_1pdt_T * ( CS%str_t(i,j) + (I_EC2 * dt_2Tdamp) * &
( zeta(i,j) * sh_Dt(i,j) ) )
enddo ; enddo
else
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,CS,I_1pdt_T,dt_2Tdamp,zeta, &
!$OMP sh_Dd,I_EC2,sh_Dt,pres_mice,mice)
do j=jsc-1,jec+1 ; do i=isc-1,iec+1
! This expression uses that Pres=2*del_sh*zeta with an elliptic yield curve.
CS%str_d(i,j) = I_1pdt_T * ( CS%str_d(i,j) + dt_2Tdamp * &
( zeta(i,j) * sh_Dd(i,j) - 0.5*pres_mice(i,j)*mice(i,j) ) )
CS%str_t(i,j) = I_1pdt_T * ( CS%str_t(i,j) + (I_EC2 * dt_2Tdamp) * &
( zeta(i,j) * sh_Dt(i,j) ) )
enddo ; enddo
endif
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,CS,I_1pdt_T,I_EC2,dt_2Tdamp, &
!$OMP G,zeta,mi_ratio_A_q,sh_Ds)
do J=jsc-1,jec ; do I=isc-1,iec
! zeta is already set to 0 over land.
CS%str_s(I,J) = I_1pdt_T * ( CS%str_s(I,J) + (I_EC2 * dt_2Tdamp) * &
( ((G%areaT(i,j)*zeta(i,j) + G%areaT(i+1,j+1)*zeta(i+1,j+1)) + &
(G%areaT(i+1,j)*zeta(i+1,j) + G%areaT(i,j+1)*zeta(i,j+1))) * &
mi_ratio_A_q(I,J) * sh_Ds(I,J) ) )
enddo ; enddo
cdRho = CS%cdw * CS%Rho_ocean
! Save the current values of u for later use in updating v.
do I=isc-1,iec
u_tmp(I,jsc-1) = ui(I,jsc-1) ; u_tmp(I,jec+1) = ui(I,jec+1) ;
enddo
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,u_tmp,ui,vi,azon,bzon,czon,dzon, &
!$OMP G,CS,dy2T,dx2B,vo,uo,Cor_u,f2dt_u,I1_f2dt2_u, &
!$OMP mi_u,dt,PFu,fxat,I_cdRhoDt,cdRho,m_neglect,fxoc, &
!$OMP fxic,fxic_d,fxic_t,fxic_s,do_trunc_its, &
!$OMP ui_min_trunc,ui_max_trunc,drag_max) &
!$OMP private(Cor,fxic_now,v2_at_u,uio_init,drag_u,b_vel0, &
!$OMP m_uio_explicit,uio_pred,uio_C)
do j=jsc,jec ; do I=isc-1,iec
! Save the current values of u for later use in updating v.
u_tmp(I,j) = ui(I,j)