-
Notifications
You must be signed in to change notification settings - Fork 0
/
MOM_regridding.F90
2472 lines (2168 loc) · 117 KB
/
MOM_regridding.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
!> Generates vertical grids as part of the ALE algorithm
module MOM_regridding
! This file is part of MOM6. See LICENSE.md for the license.
use MOM_error_handler, only : MOM_error, FATAL, WARNING
use MOM_file_parser, only : param_file_type, get_param, log_param
use MOM_io, only : file_exists, field_exists, field_size, MOM_read_data
use MOM_io, only : verify_variable_units, slasher
use MOM_unit_scaling, only : unit_scale_type
use MOM_grid, only : ocean_grid_type
use MOM_verticalGrid, only : verticalGrid_type
use MOM_EOS, only : EOS_type !, calculate_density
use MOM_string_functions, only : uppercase, extractWord, extract_integer, extract_real
use MOM_remapping, only : remapping_CS
use regrid_consts, only : state_dependent, coordinateUnits
use regrid_consts, only : coordinateMode, DEFAULT_COORDINATE_MODE
use regrid_consts, only : REGRIDDING_LAYER, REGRIDDING_ZSTAR
use regrid_consts, only : REGRIDDING_RHO, REGRIDDING_SIGMA
use regrid_consts, only : REGRIDDING_ARBITRARY, REGRIDDING_SIGMA_SHELF_ZSTAR
use regrid_consts, only : REGRIDDING_SLIGHT, REGRIDDING_ADAPTIVE, REGRIDDING_HYCOM1
use regrid_interp, only : interp_CS_type, set_interp_scheme, set_interp_extrap
use coord_zlike, only : init_coord_zlike, zlike_CS, set_zlike_params, build_zstar_column, end_coord_zlike
use coord_sigma, only : init_coord_sigma, sigma_CS, set_sigma_params, build_sigma_column, end_coord_sigma
use coord_rho, only : init_coord_rho, rho_CS, set_rho_params, build_rho_column, end_coord_rho
!use coord_rho, only : old_inflate_layers_1d
use coord_hycom, only : init_coord_hycom, hycom_CS, set_hycom_params, build_hycom1_column, end_coord_hycom
!use coord_slight, only : init_coord_slight, slight_CS, set_slight_params, build_slight_column, end_coord_slight
!use coord_adapt, only : init_coord_adapt, adapt_CS, set_adapt_params, build_adapt_column, end_coord_adapt
implicit none ; private
#include <MOM_memory.h>
! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
! vary with the Boussinesq approximation, the Boussinesq variant is given first.
type :: thermo_var_ptrs
real, dimension(:,:,:), pointer :: T=>NULL()
real, dimension(:,:,:), pointer :: S=>NULL()
type(EOS_type), pointer :: eqn_of_state => NULL()
real :: p_ref
end type thermo_var_ptrs
!> Regridding control structure
type, public :: regridding_CS
!> This array is set by function setCoordinateResolution()
!! It contains the "resolution" or delta coordinate of the target
!! coordinate. It has the units of the target coordinate, e.g.
!! [Z ~> m] for z*, non-dimensional for sigma, etc.
real, dimension(:), allocatable :: coordinateResolution
!> This is a scaling factor that restores coordinateResolution to values in
!! the natural units for output.
real :: coord_scale = 1.0
!> This array is set by function set_target_densities()
!! This array is the nominal coordinate of interfaces and is the
!! running sum of coordinateResolution, in [R ~> kg m-3]. i.e.
!! target_density(k+1) = coordinateResolution(k) + coordinateResolution(k)
!! It is only used in "rho", "SLight" or "Hycom" mode.
real, dimension(:), allocatable :: target_density
!> A flag to indicate that the target_density arrays has been filled with data.
logical :: target_density_set = .false.
!> This array is set by function set_regrid_max_depths()
!! It specifies the maximum depth that every interface is allowed to take [H ~> m or kg m-2].
real, dimension(:), allocatable :: max_interface_depths
!> This array is set by function set_regrid_max_thickness()
!! It specifies the maximum depth that every interface is allowed to take [H ~> m or kg m-2].
real, dimension(:), allocatable :: max_layer_thickness
integer :: nk !< Number of layers/levels in generated grid
!> Indicates which grid to use in the vertical (z*, sigma, target interface
!! densities)
integer :: regridding_scheme
!> Interpolation control structure
type(interp_CS_type) :: interp_CS
!> Minimum thickness allowed when building the new grid through regridding [H ~> m or kg m-2].
real :: min_thickness
!> Reference pressure for potential density calculations [R L2 T-2 ~> Pa]
real :: ref_pressure = 2.e7
!> Weight given to old coordinate when blending between new and old grids [nondim]
!! Used only below depth_of_time_filter_shallow, with a cubic variation
!! from zero to full effect between depth_of_time_filter_shallow and
!! depth_of_time_filter_deep.
real :: old_grid_weight = 0.
!> Depth above which no time-filtering of grid is applied [H ~> m or kg m-2]
real :: depth_of_time_filter_shallow = 0.
!> Depth below which time-filtering of grid is applied at full effect [H ~> m or kg m-2]
real :: depth_of_time_filter_deep = 0.
!> Fraction (between 0 and 1) of compressibility to add to potential density
!! profiles when interpolating for target grid positions. [nondim]
real :: compressibility_fraction = 0.
!> If true, each interface is given a maximum depth based on a rescaling of
!! the indexing of coordinateResolution.
logical :: set_maximum_depths = .false.
!> A scaling factor (> 1) of the rate at which the coordinateResolution list
!! is traversed to set the minimum depth of interfaces.
real :: max_depth_index_scale = 2.0
!> If true, integrate for interface positions from the top downward.
!! If false, integrate from the bottom upward, as does the rest of the model.
logical :: integrate_downward_for_e = .true.
!> If true, use the order of arithmetic and expressions that recover the remapping answers from 2018.
!! If false, use more robust forms of the same remapping expressions.
logical :: remap_answers_2018 = .true.
type(zlike_CS), pointer :: zlike_CS => null() !< Control structure for z-like coordinate generator
type(sigma_CS), pointer :: sigma_CS => null() !< Control structure for sigma coordinate generator
type(rho_CS), pointer :: rho_CS => null() !< Control structure for rho coordinate generator
type(hycom_CS), pointer :: hycom_CS => null() !< Control structure for hybrid coordinate generator
! type(slight_CS), pointer :: slight_CS => null() !< Control structure for Slight-coordinate generator
! type(adapt_CS), pointer :: adapt_CS => null() !< Control structure for adaptive coordinate generator
end type
! The following routines are visible to the outside world
public initialize_regridding, end_regridding, regridding_main, thermo_var_ptrs
!public inflate_vanished_layers_old, check_remapping_grid, check_grid_column
public check_remapping_grid, check_grid_column
public set_regrid_params, get_regrid_size
public uniformResolution, setCoordinateResolution
!public build_rho_column
public set_target_densities_from_GV, set_target_densities
public set_regrid_max_depths, set_regrid_max_thickness
public getCoordinateResolution, getCoordinateInterfaces
public getCoordinateUnits, getCoordinateShortName, getStaticThickness
public DEFAULT_COORDINATE_MODE
public get_zlike_CS, get_sigma_CS!, get_rho_CS
!> Documentation for coordinate options
character(len=*), parameter, public :: regriddingCoordinateModeDoc = &
" LAYER - Isopycnal or stacked shallow water layers\n"//&
" ZSTAR, Z* - stretched geopotential z*\n"//&
" SIGMA_SHELF_ZSTAR - stretched geopotential z* ignoring shelf\n"//&
" SIGMA - terrain following coordinates\n"//&
" RHO - continuous isopycnal\n"//&
" HYCOM1 - HyCOM-like hybrid coordinate\n"//&
" SLIGHT - stretched coordinates above continuous isopycnal\n"//&
" ADAPTIVE - optimize for smooth neutral density surfaces"
!> Documentation for regridding interpolation schemes
character(len=*), parameter, public :: regriddingInterpSchemeDoc = &
" P1M_H2 (2nd-order accurate)\n"//&
" P1M_H4 (2nd-order accurate)\n"//&
" P1M_IH4 (2nd-order accurate)\n"//&
" PLM (2nd-order accurate)\n"//&
" PPM_H4 (3rd-order accurate)\n"//&
" PPM_IH4 (3rd-order accurate)\n"//&
" P3M_IH4IH3 (4th-order accurate)\n"//&
" P3M_IH6IH5 (4th-order accurate)\n"//&
" PQM_IH4IH3 (4th-order accurate)\n"//&
" PQM_IH6IH5 (5th-order accurate)"
!> Default interpolation scheme
character(len=*), parameter, public :: regriddingDefaultInterpScheme = "P1M_H2"
!> Default mode for boundary extrapolation
logical, parameter, public :: regriddingDefaultBoundaryExtrapolation = .false.
!> Default minimum thickness for some coordinate generation modes
real, parameter, public :: regriddingDefaultMinThickness = 1.e-3
#undef __DO_SAFETY_CHECKS__
contains
!> Initialization and configures a regridding control structure based on customizable run-time parameters
subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_mode, param_prefix, param_suffix)
type(regridding_CS), intent(inout) :: CS !< Regridding control structure
type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
real, intent(in) :: max_depth !< The maximum depth of the ocean [Z ~> m].
type(param_file_type), intent(in) :: param_file !< Parameter file
character(len=*), intent(in) :: mdl !< Name of calling module.
character(len=*), intent(in) :: coord_mode !< Coordinate mode
character(len=*), intent(inout), optional :: param_prefix !< String to prefix to parameter names.
!! If empty, causes main model parameters to be used.
character(len=*), intent(inout), optional :: param_suffix !< String to append to parameter names.
! Local variables
integer :: ke ! Number of levels
character(len=80) :: string, string2, varName ! Temporary strings
character(len=40) :: coord_units, param_name, coord_res_param ! Temporary strings
character(len=200) :: inputdir, fileName
character(len=320) :: message ! Temporary strings
character(len=12) :: expected_units, alt_units ! Temporary strings
logical :: tmpLogical, fix_haloclines, set_max, do_sum, main_parameters
logical :: coord_is_state_dependent, ierr
logical :: default_2018_answers, remap_answers_2018
real :: filt_len, strat_tol, index_scale, tmpReal, P_Ref
real :: maximum_depth ! The maximum depth of the ocean [m] (not in Z).
real :: dz_fixed_sfc, Rho_avg_depth, nlay_sfc_int
real :: adaptTimeRatio, adaptZoom, adaptZoomCoeff, adaptBuoyCoeff, adaptAlpha
real :: adaptDrho0 ! Reference density difference for stratification-dependent diffusion. [R ~> kg m-3]
integer :: nz_fixed_sfc, k, nzf(4)
real, dimension(:), allocatable :: dz ! Resolution (thickness) in units of coordinate, which may be [m]
! or [Z ~> m] or [H ~> m or kg m-2] or [R ~> kg m-3] or other units.
real, dimension(:), allocatable :: h_max ! Maximum layer thicknesses [H ~> m or kg m-2]
real, dimension(:), allocatable :: z_max ! Maximum interface depths [H ~> m or kg m-2] or other
! units depending on the coordinate
real, dimension(:), allocatable :: dz_max ! Thicknesses used to find maximum interface depths
! [H ~> m or kg m-2] or other units
real, dimension(:), allocatable :: rho_target ! Target density used in HYBRID mode [kg m-3]
! Thicknesses [m] that give level centers corresponding to table 2 of WOA09
real, dimension(40) :: woa09_dz = (/ 5., 10., 10., 15., 22.5, 25., 25., 25., &
37.5, 50., 50., 75., 100., 100., 100., 100., &
100., 100., 100., 100., 100., 100., 100., 175., &
250., 375., 500., 500., 500., 500., 500., 500., &
500., 500., 500., 500., 500., 500., 500., 500. /)
call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
inputdir = slasher(inputdir)
main_parameters=.false.
if (.not.PRESENT(param_prefix)) param_prefix=''
if (.not.PRESENT(param_suffix)) param_suffix=''
if (len_trim(param_prefix)==0) main_parameters=.true.
if (main_parameters .and. len_trim(param_suffix)>0) call MOM_error(FATAL,trim(mdl)//', initialize_regridding: '// &
'Suffix provided without prefix for parameter names!')
CS%nk = 0
CS%regridding_scheme = coordinateMode(coord_mode)
coord_is_state_dependent = state_dependent(coord_mode)
maximum_depth = US%Z_to_m*max_depth
if (main_parameters) then
! Read coordinate units parameter (main model = REGRIDDING_COORDINATE_UNITS)
call get_param(param_file, mdl, "REGRIDDING_COORDINATE_UNITS", coord_units, &
"Units of the regridding coordinate.", default=coordinateUnits(coord_mode))
else
coord_units=coordinateUnits(coord_mode)
endif
if (coord_is_state_dependent) then
if (main_parameters) then
param_name = "INTERPOLATION_SCHEME"
string2 = regriddingDefaultInterpScheme
else
param_name = trim(param_prefix)//"_INTERP_SCHEME_"//trim(param_suffix)
string2 = 'PPM_H4' ! Default for diagnostics
endif
call get_param(param_file, mdl, "INTERPOLATION_SCHEME", string, &
"This sets the interpolation scheme to use to "//&
"determine the new grid. These parameters are "//&
"only relevant when REGRIDDING_COORDINATE_MODE is "//&
"set to a function of state. Otherwise, it is not "//&
"used. It can be one of the following schemes: "//&
trim(regriddingInterpSchemeDoc), default=trim(string2))
call set_regrid_params(CS, interp_scheme=string)
call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, &
"This sets the default value for the various _2018_ANSWERS parameters.", &
default=.false.)
call get_param(param_file, mdl, "REMAPPING_2018_ANSWERS", remap_answers_2018, &
"If true, use the order of arithmetic and expressions that recover the "//&
"answers from the end of 2018. Otherwise, use updated and more robust "//&
"forms of the same expressions.", default=default_2018_answers)
call set_regrid_params(CS, remap_answers_2018=remap_answers_2018)
endif
if (main_parameters .and. coord_is_state_dependent) then
call get_param(param_file, mdl, "BOUNDARY_EXTRAPOLATION", tmpLogical, &
"When defined, a proper high-order reconstruction "//&
"scheme is used within boundary cells rather "//&
"than PCM. E.g., if PPM is used for remapping, a "//&
"PPM reconstruction will also be used within "//&
"boundary cells.", default=regriddingDefaultBoundaryExtrapolation)
call set_regrid_params(CS, boundary_extrapolation=tmpLogical)
else
call set_regrid_params(CS, boundary_extrapolation=.false.)
endif
! Read coordinate configuration parameter (main model = ALE_COORDINATE_CONFIG)
if (main_parameters) then
param_name = "ALE_COORDINATE_CONFIG"
coord_res_param = "ALE_RESOLUTION"
string2 = 'UNIFORM'
else
param_name = trim(param_prefix)//"_DEF_"//trim(param_suffix)
coord_res_param = trim(param_prefix)//"_RES_"//trim(param_suffix)
string2 = 'UNIFORM'
if (maximum_depth>3000.) string2='WOA09' ! For convenience
endif
call get_param(param_file, mdl, param_name, string, &
"Determines how to specify the coordinate "//&
"resolution. Valid options are:\n"//&
" PARAM - use the vector-parameter "//trim(coord_res_param)//"\n"//&
" UNIFORM[:N] - uniformly distributed\n"//&
" FILE:string - read from a file. The string specifies\n"//&
" the filename and variable name, separated\n"//&
" by a comma or space, e.g. FILE:lev.nc,dz\n"//&
" or FILE:lev.nc,interfaces=zw\n"//&
" WOA09[:N] - the WOA09 vertical grid (approximately)\n"//&
" FNC1:string - FNC1:dz_min,H_total,power,precision\n"//&
" HYBRID:string - read from a file. The string specifies\n"//&
" the filename and two variable names, separated\n"//&
" by a comma or space, for sigma-2 and dz. e.g.\n"//&
" HYBRID:vgrid.nc,sigma2,dz",&
default=trim(string2))
message = "The distribution of vertical resolution for the target\n"//&
"grid used for Eulerian-like coordinates. For example,\n"//&
"in z-coordinate mode, the parameter is a list of level\n"//&
"thicknesses (in m). In sigma-coordinate mode, the list\n"//&
"is of non-dimensional fractions of the water column."
if (index(trim(string),'UNIFORM')==1) then
if (len_trim(string)==7) then
ke = GV%ke ! Use model nk by default
tmpReal = maximum_depth
elseif (index(trim(string),'UNIFORM:')==1 .and. len_trim(string)>8) then
! Format is "UNIFORM:N" or "UNIFORM:N,dz"
ke = extract_integer(string(9:len_trim(string)),'',1)
tmpReal = extract_real(string(9:len_trim(string)),',',2,missing_value=maximum_depth)
else
call MOM_error(FATAL,trim(mdl)//', initialize_regridding: '// &
'Unable to interpret "'//trim(string)//'".')
endif
allocate(dz(ke))
dz(:) = uniformResolution(ke, coord_mode, tmpReal, &
US%R_to_kg_m3*(GV%Rlay(1) + 0.5*(GV%Rlay(1)-GV%Rlay(min(2,ke)))), &
US%R_to_kg_m3*(GV%Rlay(ke) + 0.5*(GV%Rlay(ke)-GV%Rlay(max(ke-1,1)))) )
if (main_parameters) call log_param(param_file, mdl, "!"//coord_res_param, dz, &
trim(message), units=trim(coord_units))
elseif (trim(string)=='PARAM') then
! Read coordinate resolution (main model = ALE_RESOLUTION)
ke = GV%ke ! Use model nk by default
allocate(dz(ke))
call get_param(param_file, mdl, coord_res_param, dz, &
trim(message), units=trim(coord_units), fail_if_missing=.true.)
elseif (index(trim(string),'FILE:')==1) then
! FILE:filename,var_name is assumed to be reading level thickness variables
! FILE:filename,interfaces=var_name reads positions
if (string(6:6)=='.' .or. string(6:6)=='/') then
! If we specified "FILE:./xyz" or "FILE:/xyz" then we have a relative or absolute path
fileName = trim( extractWord(trim(string(6:80)), 1) )
else
! Otherwise assume we should look for the file in INPUTDIR
fileName = trim(inputdir) // trim( extractWord(trim(string(6:80)), 1) )
endif
if (.not. file_exists(fileName)) call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// &
"Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(string)//")")
varName = trim( extractWord(trim(string(6:)), 2) )
if (len_trim(varName)==0) then
if (field_exists(fileName,'dz')) then; varName = 'dz'
elseif (field_exists(fileName,'dsigma')) then; varName = 'dsigma'
elseif (field_exists(fileName,'ztest')) then; varName = 'ztest'
else ; call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// &
"Coordinate variable not specified and none could be guessed.")
endif
endif
! This check fails when the variable is a dimension variable! -AJA
!if (.not. field_exists(fileName,trim(varName))) call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// &
! "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")")
if (CS%regridding_scheme == REGRIDDING_SIGMA) then
expected_units = 'nondim' ; alt_units = expected_units
elseif (CS%regridding_scheme == REGRIDDING_RHO) then
expected_units = 'kg m-3' ; alt_units = expected_units
else
expected_units = 'meters' ; alt_units = 'm'
endif
if (index(trim(varName),'interfaces=')==1) then
varName=trim(varName(12:))
call verify_variable_units(filename, varName, expected_units, message, ierr, alt_units)
if (ierr) call MOM_error(FATAL, trim(mdl)//", initialize_regridding: "//&
"Unsupported format in grid definition '"//trim(filename)//"'. Error message "//trim(message))
call field_size(trim(fileName), trim(varName), nzf)
ke = nzf(1)-1
if (ke < 1) call MOM_error(FATAL, trim(mdl)//" initialize_regridding via Var "//trim(varName)//&
"in FILE "//trim(filename)//" requires at least 2 target interface values.")
if (CS%regridding_scheme == REGRIDDING_RHO) then
allocate(rho_target(ke+1))
call MOM_read_data(trim(fileName), trim(varName), rho_target)
else
allocate(dz(ke))
allocate(z_max(ke+1))
call MOM_read_data(trim(fileName), trim(varName), z_max)
dz(:) = abs(z_max(1:ke) - z_max(2:ke+1))
deallocate(z_max)
endif
else
! Assume reading resolution
call field_size(trim(fileName), trim(varName), nzf)
ke = nzf(1)
allocate(dz(ke))
call MOM_read_data(trim(fileName), trim(varName), dz)
endif
if (main_parameters .and. (ke/=GV%ke)) then
call MOM_error(FATAL,trim(mdl)//', initialize_regridding: '// &
'Mismatch in number of model levels and "'//trim(string)//'".')
endif
if (main_parameters) call log_param(param_file, mdl, "!"//coord_res_param, dz, &
trim(message), units=coordinateUnits(coord_mode))
elseif (index(trim(string),'FNC1:')==1) then
ke = GV%ke; allocate(dz(ke))
call dz_function1( trim(string(6:)), dz )
if (main_parameters) call log_param(param_file, mdl, "!"//coord_res_param, dz, &
trim(message), units=coordinateUnits(coord_mode))
elseif (index(trim(string),'RFNC1:')==1) then
! Function used for set target interface densities
ke = rho_function1( trim(string(7:)), rho_target )
elseif (index(trim(string),'HYBRID:')==1) then
ke = GV%ke; allocate(dz(ke))
! The following assumes the FILE: syntax of above but without "FILE:" in the string
allocate(rho_target(ke+1))
fileName = trim( extractWord(trim(string(8:)), 1) )
if (fileName(1:1)/='.' .and. filename(1:1)/='/') fileName = trim(inputdir) // trim( fileName )
if (.not. file_exists(fileName)) call MOM_error(FATAL,trim(mdl)//", initialize_regridding: HYBRID "// &
"Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(string)//")")
varName = trim( extractWord(trim(string(8:)), 2) )
if (.not. field_exists(fileName,varName)) call MOM_error(FATAL,trim(mdl)//", initialize_regridding: HYBRID "// &
"Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")")
call MOM_read_data(trim(fileName), trim(varName), rho_target)
varName = trim( extractWord(trim(string(8:)), 3) )
if (varName(1:5) == 'FNC1:') then ! Use FNC1 to calculate dz
call dz_function1( trim(string((index(trim(string),'FNC1:')+5):)), dz )
else ! Read dz from file
if (.not. field_exists(fileName,varName)) call MOM_error(FATAL,trim(mdl)//", initialize_regridding: HYBRID "// &
"Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")")
call MOM_read_data(trim(fileName), trim(varName), dz)
endif
if (main_parameters) then
call log_param(param_file, mdl, "!"//coord_res_param, dz, &
trim(message), units=coordinateUnits(coord_mode))
call log_param(param_file, mdl, "!TARGET_DENSITIES", rho_target, &
'HYBRID target densities for interfaces', units=coordinateUnits(coord_mode))
endif
elseif (index(trim(string),'WOA09')==1) then
if (len_trim(string)==5) then
tmpReal = 0. ; ke = 0
do while (tmpReal<maximum_depth)
ke = ke + 1
tmpReal = tmpReal + woa09_dz(ke)
enddo
elseif (index(trim(string),'WOA09:')==1) then
if (len_trim(string)==6) call MOM_error(FATAL,trim(mdl)//', initialize_regridding: '// &
'Expected string of form "WOA09:N" but got "'//trim(string)//'".')
ke = extract_integer(string(7:len_trim(string)),'',1)
endif
if (ke>40 .or. ke<1) call MOM_error(FATAL,trim(mdl)//', initialize_regridding: '// &
'For "WOA05:N" N must 0<N<41 but got "'//trim(string)//'".')
allocate(dz(ke))
dz(1:ke) = woa09_dz(1:ke)
else
call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// &
"Unrecognized coordinate configuration"//trim(string))
endif
if (main_parameters) then
! This is a work around to apparently needed to work with the from_Z initialization... ???
if (coordinateMode(coord_mode) == REGRIDDING_ZSTAR .or. &
coordinateMode(coord_mode) == REGRIDDING_HYCOM1 .or. &
coordinateMode(coord_mode) == REGRIDDING_SLIGHT .or. &
coordinateMode(coord_mode) == REGRIDDING_ADAPTIVE) then
! Adjust target grid to be consistent with maximum_depth
tmpReal = sum( dz(:) )
if (tmpReal < maximum_depth) then
dz(ke) = dz(ke) + ( maximum_depth - tmpReal )
elseif (tmpReal > maximum_depth) then
if ( dz(ke) + ( maximum_depth - tmpReal ) > 0. ) then
dz(ke) = dz(ke) + ( maximum_depth - tmpReal )
else
call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// &
"MAXIMUM_DEPTH was too shallow to adjust bottom layer of DZ!"//trim(string))
endif
endif
endif
endif
CS%nk=ke
! Target resolution (for fixed coordinates)
allocate( CS%coordinateResolution(CS%nk) ); CS%coordinateResolution(:) = -1.E30
if (state_dependent(CS%regridding_scheme)) then
! Target values
allocate( CS%target_density(CS%nk+1) ); CS%target_density(:) = -1.E30*US%kg_m3_to_R
endif
if (allocated(dz)) then
if (coordinateMode(coord_mode) == REGRIDDING_SIGMA) then
call setCoordinateResolution(dz, CS, scale=1.0)
elseif (coordinateMode(coord_mode) == REGRIDDING_RHO) then
call setCoordinateResolution(dz, CS, scale=US%kg_m3_to_R)
CS%coord_scale = US%R_to_kg_m3
elseif (coordinateMode(coord_mode) == REGRIDDING_ADAPTIVE) then
call setCoordinateResolution(dz, CS, scale=GV%m_to_H)
CS%coord_scale = GV%H_to_m
else
call setCoordinateResolution(dz, CS, scale=US%m_to_Z)
CS%coord_scale = US%Z_to_m
endif
endif
if (allocated(rho_target)) then
call set_target_densities(CS, US%kg_m3_to_R*rho_target)
print *,'Target Densities=',CS%target_density
deallocate(rho_target)
! \todo This line looks like it would overwrite the target densities set just above?
elseif (coordinateMode(coord_mode) == REGRIDDING_RHO) then
call set_target_densities_from_GV(GV, US, CS)
call log_param(param_file, mdl, "!TARGET_DENSITIES", US%R_to_kg_m3*CS%target_density(:), &
'RHO target densities for interfaces', units=coordinateUnits(coord_mode))
endif
! initialise coordinate-specific control structure
call initCoord(CS, GV, US, coord_mode)
if (main_parameters .and. coord_is_state_dependent) then
call get_param(param_file, mdl, "P_REF", P_Ref, &
"The pressure that is used for calculating the coordinate "//&
"density. (1 Pa = 1e4 dbar, so 2e7 is commonly used.) "//&
"This is only used if USE_EOS and ENABLE_THERMODYNAMICS are true.", &
units="Pa", default=2.0e7, scale=US%kg_m3_to_R*US%m_s_to_L_T**2)
call get_param(param_file, mdl, "REGRID_COMPRESSIBILITY_FRACTION", tmpReal, &
"When interpolating potential density profiles we can add "//&
"some artificial compressibility solely to make homogeneous "//&
"regions appear stratified.", units="nondim", default=0.)
call set_regrid_params(CS, compress_fraction=tmpReal, ref_pressure=P_Ref)
endif
if (main_parameters) then
call get_param(param_file, mdl, "MIN_THICKNESS", tmpReal, &
"When regridding, this is the minimum layer "//&
"thickness allowed.", units="m", scale=GV%m_to_H, &
default=regriddingDefaultMinThickness,do_not_log=.true. )
call set_regrid_params(CS, min_thickness=tmpReal)
else
call set_regrid_params(CS, min_thickness=0.)
endif
if (coordinateMode(coord_mode) == REGRIDDING_SLIGHT) then
! Set SLight-specific regridding parameters.
call get_param(param_file, mdl, "SLIGHT_DZ_SURFACE", dz_fixed_sfc, &
"The nominal thickness of fixed thickness near-surface "//&
"layers with the SLight coordinate.", units="m", default=1.0, scale=GV%m_to_H)
call get_param(param_file, mdl, "SLIGHT_NZ_SURFACE_FIXED", nz_fixed_sfc, &
"The number of fixed-depth surface layers with the SLight "//&
"coordinate.", units="nondimensional", default=2)
call get_param(param_file, mdl, "SLIGHT_SURFACE_AVG_DEPTH", Rho_avg_depth, &
"The thickness of the surface region over which to average "//&
"when calculating the density to use to define the interior "//&
"with the SLight coordinate.", units="m", default=1.0, scale=GV%m_to_H)
call get_param(param_file, mdl, "SLIGHT_NLAY_TO_INTERIOR", nlay_sfc_int, &
"The number of layers to offset the surface density when "//&
"defining where the interior ocean starts with SLight.", &
units="nondimensional", default=2.0)
call get_param(param_file, mdl, "SLIGHT_FIX_HALOCLINES", fix_haloclines, &
"If true, identify regions above the reference pressure "//&
"where the reference pressure systematically underestimates "//&
"the stratification and use this in the definition of the "//&
"interior with the SLight coordinate.", default=.false.)
call set_regrid_params(CS, dz_min_surface=dz_fixed_sfc, &
nz_fixed_surface=nz_fixed_sfc, Rho_ML_avg_depth=Rho_avg_depth, &
nlay_ML_to_interior=nlay_sfc_int, fix_haloclines=fix_haloclines)
if (fix_haloclines) then
! Set additional parameters related to SLIGHT_FIX_HALOCLINES.
call get_param(param_file, mdl, "HALOCLINE_FILTER_LENGTH", filt_len, &
"A length scale over which to smooth the temperature and "//&
"salinity before identifying erroneously unstable haloclines.", &
units="m", default=2.0, scale=GV%m_to_H)
call get_param(param_file, mdl, "HALOCLINE_STRAT_TOL", strat_tol, &
"A tolerance for the ratio of the stratification of the "//&
"apparent coordinate stratification to the actual value "//&
"that is used to identify erroneously unstable haloclines. "//&
"This ratio is 1 when they are equal, and sensible values "//&
"are between 0 and 0.5.", units="nondimensional", default=0.2)
call set_regrid_params(CS, halocline_filt_len=filt_len, &
halocline_strat_tol=strat_tol)
endif
endif
if (coordinateMode(coord_mode) == REGRIDDING_ADAPTIVE) then
call get_param(param_file, mdl, "ADAPT_TIME_RATIO", adaptTimeRatio, &
"Ratio of ALE timestep to grid timescale.", units="nondim", default=1.0e-1)
call get_param(param_file, mdl, "ADAPT_ZOOM_DEPTH", adaptZoom, &
"Depth of near-surface zooming region.", units="m", default=200.0, scale=GV%m_to_H)
call get_param(param_file, mdl, "ADAPT_ZOOM_COEFF", adaptZoomCoeff, &
"Coefficient of near-surface zooming diffusivity.", units="nondim", default=0.2)
call get_param(param_file, mdl, "ADAPT_BUOY_COEFF", adaptBuoyCoeff, &
"Coefficient of buoyancy diffusivity.", units="nondim", default=0.8)
call get_param(param_file, mdl, "ADAPT_ALPHA", adaptAlpha, &
"Scaling on optimization tendency.", units="nondim", default=1.0)
call get_param(param_file, mdl, "ADAPT_DO_MIN_DEPTH", tmpLogical, &
"If true, make a HyCOM-like mixed layer by preventing interfaces "//&
"from being shallower than the depths specified by the regridding coordinate.", &
default=.false.)
call get_param(param_file, mdl, "ADAPT_DRHO0", adaptDrho0, &
"Reference density difference for stratification-dependent diffusion.", &
units="kg m-3", default=0.5, scale=US%kg_m3_to_R)
call set_regrid_params(CS, adaptTimeRatio=adaptTimeRatio, adaptZoom=adaptZoom, &
adaptZoomCoeff=adaptZoomCoeff, adaptBuoyCoeff=adaptBuoyCoeff, adaptAlpha=adaptAlpha, &
adaptDoMin=tmpLogical, adaptDrho0=adaptDrho0)
endif
if (main_parameters .and. coord_is_state_dependent) then
call get_param(param_file, mdl, "MAXIMUM_INT_DEPTH_CONFIG", string, &
"Determines how to specify the maximum interface depths.\n"//&
"Valid options are:\n"//&
" NONE - there are no maximum interface depths\n"//&
" PARAM - use the vector-parameter MAXIMUM_INTERFACE_DEPTHS\n"//&
" FILE:string - read from a file. The string specifies\n"//&
" the filename and variable name, separated\n"//&
" by a comma or space, e.g. FILE:lev.nc,Z\n"//&
" FNC1:string - FNC1:dz_min,H_total,power,precision",&
default='NONE')
message = "The list of maximum depths for each interface."
allocate(z_max(ke+1))
allocate(dz_max(ke))
if ( trim(string) == "NONE") then
! Do nothing.
elseif ( trim(string) == "PARAM") then
call get_param(param_file, mdl, "MAXIMUM_INTERFACE_DEPTHS", z_max, &
trim(message), units="m", scale=GV%m_to_H, fail_if_missing=.true.)
call set_regrid_max_depths(CS, z_max)
elseif (index(trim(string),'FILE:')==1) then
if (string(6:6)=='.' .or. string(6:6)=='/') then
! If we specified "FILE:./xyz" or "FILE:/xyz" then we have a relative or absolute path
fileName = trim( extractWord(trim(string(6:80)), 1) )
else
! Otherwise assume we should look for the file in INPUTDIR
fileName = trim(inputdir) // trim( extractWord(trim(string(6:80)), 1) )
endif
if (.not. file_exists(fileName)) call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// &
"Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(string)//")")
do_sum = .false.
varName = trim( extractWord(trim(string(6:)), 2) )
if (.not. field_exists(fileName,varName)) call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// &
"Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")")
if (len_trim(varName)==0) then
if (field_exists(fileName,'z_max')) then; varName = 'z_max'
elseif (field_exists(fileName,'dz')) then; varName = 'dz' ; do_sum = .true.
elseif (field_exists(fileName,'dz_max')) then; varName = 'dz_max' ; do_sum = .true.
else ; call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// &
"MAXIMUM_INT_DEPTHS variable not specified and none could be guessed.")
endif
endif
if (do_sum) then
call MOM_read_data(trim(fileName), trim(varName), dz_max)
z_max(1) = 0.0 ; do K=1,ke ; z_max(K+1) = z_max(K) + dz_max(k) ; enddo
else
call MOM_read_data(trim(fileName), trim(varName), z_max)
endif
call log_param(param_file, mdl, "!MAXIMUM_INT_DEPTHS", z_max, &
trim(message), units=coordinateUnits(coord_mode))
call set_regrid_max_depths(CS, z_max, GV%m_to_H)
elseif (index(trim(string),'FNC1:')==1) then
call dz_function1( trim(string(6:)), dz_max )
if ((coordinateMode(coord_mode) == REGRIDDING_SLIGHT) .and. &
(dz_fixed_sfc > 0.0)) then
do k=1,nz_fixed_sfc ; dz_max(k) = dz_fixed_sfc ; enddo
endif
z_max(1) = 0.0 ; do K=1,ke ; z_max(K+1) = z_max(K) + dz_max(K) ; enddo
call log_param(param_file, mdl, "!MAXIMUM_INT_DEPTHS", z_max, &
trim(message), units=coordinateUnits(coord_mode))
call set_regrid_max_depths(CS, z_max, GV%m_to_H)
else
call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// &
"Unrecognized MAXIMUM_INT_DEPTH_CONFIG "//trim(string))
endif
deallocate(z_max)
deallocate(dz_max)
! Optionally specify maximum thicknesses for each layer, enforced by moving
! the interface below a layer downward.
call get_param(param_file, mdl, "MAX_LAYER_THICKNESS_CONFIG", string, &
"Determines how to specify the maximum layer thicknesses.\n"//&
"Valid options are:\n"//&
" NONE - there are no maximum layer thicknesses\n"//&
" PARAM - use the vector-parameter MAX_LAYER_THICKNESS\n"//&
" FILE:string - read from a file. The string specifies\n"//&
" the filename and variable name, separated\n"//&
" by a comma or space, e.g. FILE:lev.nc,Z\n"//&
" FNC1:string - FNC1:dz_min,H_total,power,precision",&
default='NONE')
message = "The list of maximum thickness for each layer."
allocate(h_max(ke))
if ( trim(string) == "NONE") then
! Do nothing.
elseif ( trim(string) == "PARAM") then
call get_param(param_file, mdl, "MAX_LAYER_THICKNESS", h_max, &
trim(message), units="m", fail_if_missing=.true., scale=GV%m_to_H)
call set_regrid_max_thickness(CS, h_max)
elseif (index(trim(string),'FILE:')==1) then
if (string(6:6)=='.' .or. string(6:6)=='/') then
! If we specified "FILE:./xyz" or "FILE:/xyz" then we have a relative or absolute path
fileName = trim( extractWord(trim(string(6:80)), 1) )
else
! Otherwise assume we should look for the file in INPUTDIR
fileName = trim(inputdir) // trim( extractWord(trim(string(6:80)), 1) )
endif
if (.not. file_exists(fileName)) call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// &
"Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(string)//")")
varName = trim( extractWord(trim(string(6:)), 2) )
if (.not. field_exists(fileName,varName)) call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// &
"Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")")
if (len_trim(varName)==0) then
if (field_exists(fileName,'h_max')) then; varName = 'h_max'
elseif (field_exists(fileName,'dz_max')) then; varName = 'dz_max'
else ; call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// &
"MAXIMUM_INT_DEPTHS variable not specified and none could be guessed.")
endif
endif
call MOM_read_data(trim(fileName), trim(varName), h_max)
call log_param(param_file, mdl, "!MAX_LAYER_THICKNESS", h_max, &
trim(message), units=coordinateUnits(coord_mode))
call set_regrid_max_thickness(CS, h_max, GV%m_to_H)
elseif (index(trim(string),'FNC1:')==1) then
call dz_function1( trim(string(6:)), h_max )
call log_param(param_file, mdl, "!MAX_LAYER_THICKNESS", h_max, &
trim(message), units=coordinateUnits(coord_mode))
call set_regrid_max_thickness(CS, h_max, GV%m_to_H)
else
call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// &
"Unrecognized MAX_LAYER_THICKNESS_CONFIG "//trim(string))
endif
deallocate(h_max)
endif
if (allocated(dz)) deallocate(dz)
end subroutine initialize_regridding
!> Deallocation of regridding memory
subroutine end_regridding(CS)
type(regridding_CS), intent(inout) :: CS !< Regridding control structure
if (associated(CS%zlike_CS)) call end_coord_zlike(CS%zlike_CS)
if (associated(CS%sigma_CS)) call end_coord_sigma(CS%sigma_CS)
! if (associated(CS%rho_CS)) call end_coord_rho(CS%rho_CS)
if (associated(CS%hycom_CS)) call end_coord_hycom(CS%hycom_CS)
! if (associated(CS%slight_CS)) call end_coord_slight(CS%slight_CS)
! if (associated(CS%adapt_CS)) call end_coord_adapt(CS%adapt_CS)
deallocate( CS%coordinateResolution )
if (allocated(CS%target_density)) deallocate( CS%target_density )
if (allocated(CS%max_interface_depths) ) deallocate( CS%max_interface_depths )
if (allocated(CS%max_layer_thickness) ) deallocate( CS%max_layer_thickness )
end subroutine end_regridding
!------------------------------------------------------------------------------
!> Dispatching regridding routine for orchestrating regridding & remapping
subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface, frac_shelf_h, conv_adjust)
!------------------------------------------------------------------------------
! This routine takes care of (1) building a new grid and (2) remapping between
! the old grid and the new grid. The creation of the new grid can be based
! on z coordinates, target interface densities, sigma coordinates or any
! arbitrary coordinate system.
! The MOM6 interface positions are always calculated from the bottom up by
! accumulating the layer thicknesses starting at z=-G%bathyT. z increases
! upwards (decreasing k-index).
! The new grid is defined by the change in position of those interfaces in z
! dzInterface = zNew - zOld.
! Thus, if the regridding inflates the top layer, hNew(1) > hOld(1), then the
! second interface moves downward, zNew(2) < zOld(2), and dzInterface(2) < 0.
! hNew(k) = hOld(k) - dzInterface(k+1) + dzInterface(k)
! IMPORTANT NOTE:
! This is the converse of the sign convention used in the remapping code!
!------------------------------------------------------------------------------
! Arguments
type(remapping_CS), intent(in) :: remapCS !< Remapping parameters and options
type(regridding_CS), intent(in) :: CS !< Regridding control structure
type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
real, dimension(SZI_(G),SZJ_(G), SZK_(GV)), intent(inout) :: h !< Current 3D grid obtained after
!! the last time step
type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamical variables (T, S, ...)
real, dimension(SZI_(G),SZJ_(G), CS%nk), intent(inout) :: h_new !< New 3D grid consistent with target coordinate
real, dimension(SZI_(G),SZJ_(G), CS%nk+1), intent(inout) :: dzInterface !< The change in position of each interface
real, dimension(SZI_(G),SZJ_(G)), optional, intent(in ) :: frac_shelf_h !< Fractional ice shelf coverage
logical, optional, intent(in ) :: conv_adjust !< If true, do convective adjustment
! Local variables
real :: trickGnuCompiler
logical :: use_ice_shelf
logical :: do_convective_adjustment
do_convective_adjustment = .true.
if (present(conv_adjust)) do_convective_adjustment = conv_adjust
use_ice_shelf = present(frac_shelf_h)
select case ( CS%regridding_scheme )
case ( REGRIDDING_ZSTAR )
call build_zstar_grid( CS, G, GV, h, dzInterface, frac_shelf_h )
call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new)
case ( REGRIDDING_SIGMA_SHELF_ZSTAR)
call build_zstar_grid( CS, G, GV, h, dzInterface, frac_shelf_h )
call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new)
case ( REGRIDDING_SIGMA )
call build_sigma_grid( CS, G, GV, h, dzInterface )
call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new)
case ( REGRIDDING_RHO )
if (do_convective_adjustment) call convective_adjustment(G, GV, h, tv)
call build_rho_grid( G, GV, h, tv, dzInterface, remapCS, CS, frac_shelf_h )
call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new)
! case ( REGRIDDING_ARBITRARY )
! call build_grid_arbitrary( G, GV, h, dzInterface, trickGnuCompiler, CS )
! call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new)
case ( REGRIDDING_HYCOM1 )
call build_grid_HyCOM1( G, GV, G%US, h, tv, h_new, dzInterface, CS, frac_shelf_h )
! case ( REGRIDDING_SLIGHT )
! call build_grid_SLight( G, GV, G%US, h, tv, dzInterface, CS )
! call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new)
! case ( REGRIDDING_ADAPTIVE )
! call build_grid_adaptive(G, GV, G%US, h, tv, dzInterface, remapCS, CS)
! call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new)
case default
call MOM_error(FATAL,'MOM_regridding, regridding_main: '//&
'Unknown regridding scheme selected!')
end select ! type of grid
#ifdef __DO_SAFETY_CHECKS__
call check_remapping_grid(G, GV, h, dzInterface,'in regridding_main')
#endif
end subroutine regridding_main
!> Calculates h_new from h + delta_k dzInterface
subroutine calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new)
type(regridding_CS), intent(in) :: CS !< Regridding control structure
type(ocean_grid_type), intent(in) :: G !< Grid structure
type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Old layer thicknesses (arbitrary units)
real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(in) :: dzInterface !< Change in interface positions (same as h)
real, dimension(SZI_(G),SZJ_(G),CS%nk), intent(inout) :: h_new !< New layer thicknesses (same as h)
! Local variables
integer :: i, j, k, nki
nki = min(CS%nk, GV%ke)
!$OMP parallel do default(shared)
do j = G%jsc,G%jec
do i = G%isc,G%iec
! if (G%mask2dT(i,j)>0.) then
do k=1,nki
h_new(i,j,k) = max( 0., h(i,j,k) + ( dzInterface(i,j,k) - dzInterface(i,j,k+1) ) )
enddo
if (CS%nk > GV%ke) then
do k=nki+1, CS%nk
h_new(i,j,k) = max( 0., dzInterface(i,j,k) - dzInterface(i,j,k+1) )
enddo
endif
! else
! h_new(i,j,1:nki) = h(i,j,1:nki)
! if (CS%nk > GV%ke) h_new(i,j,nki+1:CS%nk) = 0.
! ! On land points, why are we keeping the original h rather than setting to zero? -AJA
! endif
enddo
enddo
end subroutine calc_h_new_by_dz
!> Check that the total thickness of two grids match
subroutine check_remapping_grid( G, GV, h, dzInterface, msg )
type(ocean_grid_type), intent(in) :: G !< Grid structure
type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: dzInterface !< Change in interface positions
!! [H ~> m or kg m-2]
character(len=*), intent(in) :: msg !< Message to append to errors
! Local variables
integer :: i, j
!$OMP parallel do default(shared)
do j = G%jsc,G%jec ; do i = G%isc,G%iec
!if (G%mask2dT(i,j)>0.)
call check_grid_column( GV%ke, GV%Z_to_H*G%bathyT(i,j), h(i,j,:), dzInterface(i,j,:), msg )
enddo ; enddo
end subroutine check_remapping_grid
!> Check that the total thickness of new and old grids are consistent
subroutine check_grid_column( nk, depth, h, dzInterface, msg )
integer, intent(in) :: nk !< Number of cells
real, intent(in) :: depth !< Depth of bottom [Z ~> m] or arbitrary units
real, dimension(nk), intent(in) :: h !< Cell thicknesses [Z ~> m] or arbitrary units
real, dimension(nk+1), intent(in) :: dzInterface !< Change in interface positions (same units as h)
character(len=*), intent(in) :: msg !< Message to append to errors
! Local variables
integer :: k
real :: eps, total_h_old, total_h_new, h_new, z_old, z_new
eps =1. ; eps = epsilon(eps)
! Total thickness of grid h
total_h_old = 0.
do k = 1,nk
total_h_old = total_h_old + h(k)
enddo
! Integrate upwards for the interfaces consistent with the rest of MOM6
z_old = - depth
if (depth == 0.) z_old = - total_h_old
total_h_new = 0.
do k = nk,1,-1
z_old = z_old + h(k) ! Old interface position above layer k
z_new = z_old + dzInterface(k) ! New interface position based on dzInterface
h_new = h(k) + ( dzInterface(k) - dzInterface(k+1) ) ! New thickness
if (h_new<0.) then
write(0,*) 'k,h,hnew=',k,h(k),h_new
write(0,*) 'dzI(k+1),dzI(k)=',dzInterface(k+1),dzInterface(k)
call MOM_error( FATAL, 'MOM_regridding, check_grid_column: '//&
'Negative layer thickness implied by re-gridding, '//trim(msg))
endif
total_h_new = total_h_new + h_new
enddo
! Conservation by implied h_new
if (abs(total_h_new-total_h_old)>real(nk-1)*0.5*(total_h_old+total_h_new)*eps) then
write(0,*) 'nk=',nk
do k = 1,nk
write(0,*) 'k,h,hnew=',k,h(k),h(k)+(dzInterface(k)-dzInterface(k+1))
enddo
write(0,*) 'Hold,Hnew,Hnew-Hold=',total_h_old,total_h_new,total_h_new-total_h_old
write(0,*) 'eps,(n)/2*eps*H=',eps,real(nk-1)*0.5*(total_h_old+total_h_new)*eps
call MOM_error( FATAL, 'MOM_regridding, check_grid_column: '//&
'Re-gridding did NOT conserve total thickness to within roundoff '//trim(msg))
endif
! Check that the top and bottom are intentionally moving
if (dzInterface(1) /= 0.) call MOM_error( FATAL, &
'MOM_regridding, check_grid_column: Non-zero dzInterface at surface! '//trim(msg))
if (dzInterface(nk+1) /= 0.) call MOM_error( FATAL, &
'MOM_regridding, check_grid_column: Non-zero dzInterface at bottom! '//trim(msg))
end subroutine check_grid_column
!> Returns the change in interface position motion after filtering and
!! assuming the top and bottom interfaces do not move. The filtering is
!! a function of depth, and is applied as the integrated average filtering
!! over the trajectory of the interface. By design, this code can not give
!! tangled interfaces provided that z_old and z_new are not already tangled.
subroutine filtered_grid_motion( CS, nk, z_old, z_new, dz_g )
type(regridding_CS), intent(in) :: CS !< Regridding control structure
integer, intent(in) :: nk !< Number of cells in source grid
real, dimension(nk+1), intent(in) :: z_old !< Old grid position [H ~> m or kg m-2]
real, dimension(CS%nk+1), intent(in) :: z_new !< New grid position [H ~> m or kg m-2]
real, dimension(CS%nk+1), intent(inout) :: dz_g !< Change in interface positions [H ~> m or kg m-2]
! Local variables
real :: sgn ! The sign convention for downward.
real :: dz_tgt, zr1, z_old_k
real :: Aq, Bq, dz0, z0, F0
real :: zs, zd, dzwt, Idzwt
real :: wtd, Iwtd
real :: Int_zs, Int_zd, dInt_zs_zd
! For debugging:
real, dimension(nk+1) :: z_act
! real, dimension(nk+1) :: ddz_g_s, ddz_g_d
logical :: debug = .false.
integer :: k
if ((z_old(nk+1) - z_old(1)) * (z_new(CS%nk+1) - z_new(1)) < 0.0) then
call MOM_error(FATAL, "filtered_grid_motion: z_old and z_new use different sign conventions.")
elseif ((z_old(nk+1) - z_old(1)) * (z_new(CS%nk+1) - z_new(1)) == 0.0) then
! This is a massless column, so do nothing and return.
do k=1,CS%nk+1 ; dz_g(k) = 0.0 ; enddo ; return
elseif ((z_old(nk+1) - z_old(1)) + (z_new(CS%nk+1) - z_new(1)) > 0.0) then
sgn = 1.0
else
sgn = -1.0
endif
if (debug) then
do k=2,CS%nk+1
if (sgn*(z_new(k)-z_new(k-1)) < -5e-16*(abs(z_new(k))+abs(z_new(k-1))) ) &
call MOM_error(FATAL, "filtered_grid_motion: z_new is tangled.")
enddo
do k=2,nk+1
if (sgn*(z_old(k)-z_old(k-1)) < -5e-16*(abs(z_old(k))+abs(z_old(k-1))) ) &
call MOM_error(FATAL, "filtered_grid_motion: z_old is tangled.")
enddo
! ddz_g_s(:) = 0.0 ; ddz_g_d(:) = 0.0
endif
zs = CS%depth_of_time_filter_shallow
zd = CS%depth_of_time_filter_deep
wtd = 1.0 - CS%old_grid_weight
Iwtd = 1.0 / wtd
dzwt = (zd - zs)
Idzwt = 0.0 ; if (abs(zd - zs) > 0.0) Idzwt = 1.0 / (zd - zs)
dInt_zs_zd = 0.5*(1.0 + Iwtd) * (zd - zs)
Aq = 0.5*(Iwtd - 1.0)
dz_g(1) = 0.0