-
Notifications
You must be signed in to change notification settings - Fork 209
/
MOM_vert_friction.F90
2043 lines (1831 loc) · 101 KB
/
MOM_vert_friction.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
!> Implements vertical viscosity (vertvisc)
module MOM_vert_friction
! This file is part of MOM6. See LICENSE.md for the license.
use MOM_domains, only : pass_var, To_All, Omit_corners
use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
use MOM_diag_mediator, only : post_product_u, post_product_sum_u
use MOM_diag_mediator, only : post_product_v, post_product_sum_v
use MOM_diag_mediator, only : diag_ctrl
use MOM_debugging, only : uvchksum, hchksum
use MOM_error_handler, only : MOM_error, FATAL, WARNING, NOTE
use MOM_file_parser, only : get_param, log_version, param_file_type
use MOM_forcing_type, only : mech_forcing
use MOM_get_input, only : directories
use MOM_grid, only : ocean_grid_type
use MOM_open_boundary, only : ocean_OBC_type, OBC_SIMPLE, OBC_NONE, OBC_DIRECTION_E
use MOM_open_boundary, only : OBC_DIRECTION_W, OBC_DIRECTION_N, OBC_DIRECTION_S
use MOM_PointAccel, only : write_u_accel, write_v_accel, PointAccel_init
use MOM_PointAccel, only : PointAccel_CS
use MOM_time_manager, only : time_type, time_type_to_real, operator(-)
use MOM_unit_scaling, only : unit_scale_type
use MOM_variables, only : thermo_var_ptrs, vertvisc_type
use MOM_variables, only : cont_diag_ptrs, accel_diag_ptrs
use MOM_variables, only : ocean_internal_state
use MOM_verticalGrid, only : verticalGrid_type
use MOM_wave_interface, only : wave_parameters_CS
implicit none ; private
#include <MOM_memory.h>
public vertvisc, vertvisc_remnant, vertvisc_coef
public vertvisc_limit_vel, vertvisc_init, vertvisc_end
public updateCFLtruncationValue
! 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.
!> The control structure with parameters and memory for the MOM_vert_friction module
type, public :: vertvisc_CS ; private
logical :: initialized = .false. !< True if this control structure has been initialized.
real :: Hmix !< The mixed layer thickness in thickness units [H ~> m or kg m-2].
real :: Hmix_stress !< The mixed layer thickness over which the wind
!! stress is applied with direct_stress [H ~> m or kg m-2].
real :: Kvml !< The mixed layer vertical viscosity [Z2 T-1 ~> m2 s-1].
real :: Kv !< The interior vertical viscosity [Z2 T-1 ~> m2 s-1].
real :: Hbbl !< The static bottom boundary layer thickness [H ~> m or kg m-2].
real :: Kvbbl !< The vertical viscosity in the bottom boundary
!! layer [Z2 T-1 ~> m2 s-1].
real :: maxvel !< Velocity components greater than maxvel are truncated [L T-1 ~> m s-1].
real :: vel_underflow !< Velocity components smaller than vel_underflow
!! are set to 0 [L T-1 ~> m s-1].
logical :: CFL_based_trunc !< If true, base truncations on CFL numbers, not
!! absolute velocities.
real :: CFL_trunc !< Velocity components will be truncated when they
!! are large enough that the corresponding CFL number
!! exceeds this value [nondim].
real :: CFL_report !< The value of the CFL number that will cause the
!! accelerations to be reported [nondim]. CFL_report
!! will often equal CFL_trunc.
real :: truncRampTime !< The time-scale over which to ramp up the value of
!! CFL_trunc from CFL_truncS to CFL_truncE [T ~> s]
real :: CFL_truncS !< The start value of CFL_trunc
real :: CFL_truncE !< The end/target value of CFL_trunc
logical :: CFLrampingIsActivated = .false. !< True if the ramping has been initialized
type(time_type) :: rampStartTime !< The time at which the ramping of CFL_trunc starts
real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_,NK_INTERFACE_) :: &
a_u !< The u-drag coefficient across an interface [Z T-1 ~> m s-1].
real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: &
h_u !< The effective layer thickness at u-points [H ~> m or kg m-2].
real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_,NK_INTERFACE_) :: &
a_v !< The v-drag coefficient across an interface [Z T-1 ~> m s-1].
real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: &
h_v !< The effective layer thickness at v-points [H ~> m or kg m-2].
real, pointer, dimension(:,:) :: a1_shelf_u => NULL() !< The u-momentum coupling coefficient under
!! ice shelves [Z T-1 ~> m s-1]. Retained to determine stress under shelves.
real, pointer, dimension(:,:) :: a1_shelf_v => NULL() !< The v-momentum coupling coefficient under
!! ice shelves [Z T-1 ~> m s-1]. Retained to determine stress under shelves.
logical :: split !< If true, use the split time stepping scheme.
logical :: bottomdraglaw !< If true, the bottom stress is calculated with a
!! drag law c_drag*|u|*u. The velocity magnitude
!! may be an assumed value or it may be based on the
!! actual velocity in the bottommost HBBL, depending
!! on whether linear_drag is true.
logical :: Channel_drag !< If true, the drag is exerted directly on each
!! layer according to what fraction of the bottom
!! they overlie.
logical :: harmonic_visc !< If true, the harmonic mean thicknesses are used
!! to calculate the viscous coupling between layers
!! except near the bottom. Otherwise the arithmetic
!! mean thickness is used except near the bottom.
real :: harm_BL_val !< A scale to determine when water is in the boundary
!! layers based solely on harmonic mean thicknesses
!! for the purpose of determining the extent to which
!! the thicknesses used in the viscosities are upwinded.
logical :: direct_stress !< If true, the wind stress is distributed over the
!! topmost Hmix_stress of fluid and KVML may be very small.
logical :: dynamic_viscous_ML !< If true, use the results from a dynamic
!! calculation, perhaps based on a bulk Richardson
!! number criterion, to determine the mixed layer
!! thickness for viscosity.
logical :: answers_2018 !< If true, use the order of arithmetic and expressions that recover the
!! answers from the end of 2018. Otherwise, use expressions that do not
!! use an arbitrary and hard-coded maximum viscous coupling coefficient
!! between layers.
logical :: debug !< If true, write verbose checksums for debugging purposes.
integer :: nkml !< The number of layers in the mixed layer.
integer, pointer :: ntrunc !< The number of times the velocity has been
!! truncated since the last call to write_energy.
character(len=200) :: u_trunc_file !< The complete path to a file in which a column of
!! u-accelerations are written if velocity truncations occur.
character(len=200) :: v_trunc_file !< The complete path to a file in which a column of
!! v-accelerations are written if velocity truncations occur.
logical :: StokesMixing !< If true, do Stokes drift mixing via the Lagrangian current
!! (Eulerian plus Stokes drift). False by default and set
!! via STOKES_MIXING_COMBINED.
type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the
!! timing of diagnostic output.
!>@{ Diagnostic identifiers
integer :: id_du_dt_visc = -1, id_dv_dt_visc = -1, id_au_vv = -1, id_av_vv = -1
integer :: id_du_dt_str = -1, id_dv_dt_str = -1
integer :: id_h_u = -1, id_h_v = -1, id_hML_u = -1 , id_hML_v = -1
integer :: id_taux_bot = -1, id_tauy_bot = -1
integer :: id_Kv_slow = -1, id_Kv_u = -1, id_Kv_v = -1
! integer :: id_hf_du_dt_visc = -1, id_hf_dv_dt_visc = -1
integer :: id_h_du_dt_visc = -1, id_h_dv_dt_visc = -1
integer :: id_hf_du_dt_visc_2d = -1, id_hf_dv_dt_visc_2d = -1
integer :: id_h_du_dt_str = -1, id_h_dv_dt_str = -1
integer :: id_du_dt_str_visc_rem = -1, id_dv_dt_str_visc_rem = -1
!>@}
type(PointAccel_CS), pointer :: PointAccel_CSp => NULL() !< A pointer to the control structure
!! for recording accelerations leading to velocity truncations
end type vertvisc_CS
contains
!> Perform a fully implicit vertical diffusion
!! of momentum. Stress top and bottom boundary conditions are used.
!!
!! This is solving the tridiagonal system
!! \f[ \left(h_k + a_{k + 1/2} + a_{k - 1/2} + r_k\right) u_k^{n+1}
!! = h_k u_k^n + a_{k + 1/2} u_{k+1}^{n+1} + a_{k - 1/2} u_{k-1}^{n+1} \f]
!! where \f$a_{k + 1/2} = \Delta t \nu_{k + 1/2} / h_{k + 1/2}\f$
!! is the <em>interfacial coupling thickness per time step</em>,
!! encompassing background viscosity as well as contributions from
!! enhanced mixed and bottom layer viscosities.
!! $r_k$ is a Rayleight drag term due to channel drag.
!! There is an additional stress term on the right-hand side
!! if DIRECT_STRESS is true, applied to the surface layer.
subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, &
taux_bot, tauy_bot, Waves)
type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
intent(inout) :: u !< Zonal velocity [L T-1 ~> m s-1]
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
intent(inout) :: v !< Meridional velocity [L T-1 ~> m s-1]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
type(vertvisc_type), intent(inout) :: visc !< Viscosities and bottom drag
real, intent(in) :: dt !< Time increment [T ~> s]
type(ocean_OBC_type), pointer :: OBC !< Open boundary condition structure
type(accel_diag_ptrs), intent(inout) :: ADp !< Accelerations in the momentum
!! equations for diagnostics
type(cont_diag_ptrs), intent(inout) :: CDp !< Continuity equation terms
type(vertvisc_CS), pointer :: CS !< Vertical viscosity control structure
real, dimension(SZIB_(G),SZJ_(G)), &
optional, intent(out) :: taux_bot !< Zonal bottom stress from ocean to
!! rock [R L Z T-2 ~> Pa]
real, dimension(SZI_(G),SZJB_(G)), &
optional, intent(out) :: tauy_bot !< Meridional bottom stress from ocean to
!! rock [R L Z T-2 ~> Pa]
type(wave_parameters_CS), &
optional, pointer :: Waves !< Container for wave/Stokes information
! Fields from forces used in this subroutine:
! taux: Zonal wind stress [R L Z T-2 ~> Pa].
! tauy: Meridional wind stress [R L Z T-2 ~> Pa].
! Local variables
real :: b1(SZIB_(G)) ! A variable used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1].
real :: c1(SZIB_(G),SZK_(GV)) ! A variable used by the tridiagonal solver [nondim].
real :: d1(SZIB_(G)) ! d1=1-c1 is used by the tridiagonal solver [nondim].
real :: Ray(SZIB_(G),SZK_(GV)) ! Ray is the Rayleigh-drag velocity [Z T-1 ~> m s-1].
real :: b_denom_1 ! The first term in the denominator of b1 [H ~> m or kg m-2].
real :: Hmix ! The mixed layer thickness over which stress
! is applied with direct_stress [H ~> m or kg m-2].
real :: I_Hmix ! The inverse of Hmix [H-1 ~> m-1 or m2 kg-1].
real :: Idt ! The inverse of the time step [T-1 ~> s-1].
real :: dt_Rho0 ! The time step divided by the mean density [T H Z-1 R-1 ~> s m3 kg-1 or s].
real :: dt_Z_to_H ! The time step times the conversion from Z to the
! units of thickness - [T H Z-1 ~> s or s kg m-3].
real :: h_neglect ! A thickness that is so small it is usually lost
! in roundoff and can be neglected [H ~> m or kg m-2].
real :: stress ! The surface stress times the time step, divided
! by the density [H L T-1 ~> m2 s-1 or kg m-1 s-1].
real :: accel_underflow ! An acceleration magnitude that is so small that values that are less
! than this are diagnosed as 0 [L T-2 ~> m s-2].
real :: zDS, hfr, h_a ! Temporary variables used with direct_stress.
real :: surface_stress(SZIB_(G))! The same as stress, unless the wind stress
! stress is applied as a body force [H L T-1 ~> m2 s-1 or kg m-1 s-1].
logical :: do_i(SZIB_(G))
logical :: DoStokesMixing
integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, n
is = G%isc ; ie = G%iec; js = G%jsc; je = G%jec
Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB ; nz = GV%ke
if (.not.associated(CS)) call MOM_error(FATAL,"MOM_vert_friction(visc): "// &
"Module must be initialized before it is used.")
if (.not.CS%initialized) call MOM_error(FATAL,"MOM_vert_friction(visc): "// &
"Module must be initialized before it is used.")
if (CS%direct_stress) then
Hmix = CS%Hmix_stress
I_Hmix = 1.0 / Hmix
endif
dt_Rho0 = dt / GV%H_to_RZ
dt_Z_to_H = dt*GV%Z_to_H
h_neglect = GV%H_subroundoff
Idt = 1.0 / dt
accel_underflow = CS%vel_underflow * Idt
!Check if Stokes mixing allowed if requested (present and associated)
DoStokesMixing=.false.
if (CS%StokesMixing) then
if (present(Waves)) DoStokesMixing = associated(Waves)
if (.not. DoStokesMixing) &
call MOM_error(FATAL,"Stokes Mixing called without allocated"//&
"Waves Control Structure")
endif
do k=1,nz ; do i=Isq,Ieq ; Ray(i,k) = 0.0 ; enddo ; enddo
! Update the zonal velocity component using a modification of a standard
! tridagonal solver.
!$OMP parallel do default(shared) firstprivate(Ray) &
!$OMP private(do_i,surface_stress,zDS,stress,h_a,hfr, &
!$OMP b_denom_1,b1,d1,c1)
do j=G%jsc,G%jec
do I=Isq,Ieq ; do_i(I) = (G%mask2dCu(I,j) > 0) ; enddo
! When mixing down Eulerian current + Stokes drift add before calling solver
if (DoStokesMixing) then ; do k=1,nz ; do I=Isq,Ieq
if (do_i(I)) u(I,j,k) = u(I,j,k) + Waves%Us_x(I,j,k)
enddo ; enddo ; endif
if (associated(ADp%du_dt_visc)) then ; do k=1,nz ; do I=Isq,Ieq
ADp%du_dt_visc(I,j,k) = u(I,j,k)
enddo ; enddo ; endif
if (associated(ADp%du_dt_str)) then ; do k=1,nz ; do I=Isq,Ieq
ADp%du_dt_str(I,j,k) = 0.0
enddo ; enddo ; endif
! One option is to have the wind stress applied as a body force
! over the topmost Hmix fluid. If DIRECT_STRESS is not defined,
! the wind stress is applied as a stress boundary condition.
if (CS%direct_stress) then
do I=Isq,Ieq ; if (do_i(I)) then
surface_stress(I) = 0.0
zDS = 0.0
stress = dt_Rho0 * forces%taux(I,j)
do k=1,nz
h_a = 0.5 * (h(I,j,k) + h(I+1,j,k)) + h_neglect
hfr = 1.0 ; if ((zDS+h_a) > Hmix) hfr = (Hmix - zDS) / h_a
u(I,j,k) = u(I,j,k) + I_Hmix * hfr * stress
if (associated(ADp%du_dt_str)) ADp%du_dt_str(i,J,k) = (I_Hmix * hfr * stress) * Idt
zDS = zDS + h_a ; if (zDS >= Hmix) exit
enddo
endif ; enddo ! end of i loop
else ; do I=Isq,Ieq
surface_stress(I) = dt_Rho0 * (G%mask2dCu(I,j)*forces%taux(I,j))
enddo ; endif ! direct_stress
if (CS%Channel_drag) then ; do k=1,nz ; do I=Isq,Ieq
Ray(I,k) = visc%Ray_u(I,j,k)
enddo ; enddo ; endif
! perform forward elimination on the tridiagonal system
!
! denote the diagonal of the system as b_k, the subdiagonal as a_k
! and the superdiagonal as c_k. The right-hand side terms are d_k.
!
! ignoring the rayleigh drag contribution,
! we have a_k = -dt_Z_to_H * a_u(k)
! b_k = h_u(k) + dt_Z_to_H * (a_u(k) + a_u(k+1))
! c_k = -dt_Z_to_H * a_u(k+1)
!
! for forward elimination, we want to:
! calculate c'_k = - c_k / (b_k + a_k c'_(k-1))
! and d'_k = (d_k - a_k d'_(k-1)) / (b_k + a_k c'_(k-1))
! where c'_1 = c_1/b_1 and d'_1 = d_1/b_1
! (see Thomas' tridiagonal matrix algorithm)
!
! b1 is the denominator term 1 / (b_k + a_k c'_(k-1))
! b_denom_1 is (b_k + a_k + c_k) - a_k(1 - c'_(k-1))
! = (b_k + c_k + c'_(k-1))
! this is done so that d1 = b1 * b_denom_1 = 1 - c'_(k-1)
! c1(k) is -c'_(k - 1)
! and the right-hand-side is destructively updated to be d'_k
!
do I=Isq,Ieq ; if (do_i(I)) then
b_denom_1 = CS%h_u(I,j,1) + dt_Z_to_H * (Ray(I,1) + CS%a_u(I,j,1))
b1(I) = 1.0 / (b_denom_1 + dt_Z_to_H*CS%a_u(I,j,2))
d1(I) = b_denom_1 * b1(I)
u(I,j,1) = b1(I) * (CS%h_u(I,j,1) * u(I,j,1) + surface_stress(I))
if (associated(ADp%du_dt_str)) &
ADp%du_dt_str(I,j,1) = b1(I) * (CS%h_u(I,j,1) * ADp%du_dt_str(I,j,1) + surface_stress(I)*Idt)
endif ; enddo
do k=2,nz ; do I=Isq,Ieq ; if (do_i(I)) then
c1(I,k) = dt_Z_to_H * CS%a_u(I,j,K) * b1(I)
b_denom_1 = CS%h_u(I,j,k) + dt_Z_to_H * (Ray(I,k) + CS%a_u(I,j,K)*d1(I))
b1(I) = 1.0 / (b_denom_1 + dt_Z_to_H * CS%a_u(I,j,K+1))
d1(I) = b_denom_1 * b1(I)
u(I,j,k) = (CS%h_u(I,j,k) * u(I,j,k) + &
dt_Z_to_H * CS%a_u(I,j,K) * u(I,j,k-1)) * b1(I)
if (associated(ADp%du_dt_str)) &
ADp%du_dt_str(I,j,k) = (CS%h_u(I,j,k) * ADp%du_dt_str(I,j,k) + &
dt_Z_to_H * CS%a_u(I,j,K) * ADp%du_dt_str(I,j,k-1)) * b1(I)
endif ; enddo ; enddo
! back substitute to solve for the new velocities
! u_k = d'_k - c'_k x_(k+1)
do k=nz-1,1,-1 ; do I=Isq,Ieq ; if (do_i(I)) then
u(I,j,k) = u(I,j,k) + c1(I,k+1) * u(I,j,k+1)
endif ; enddo ; enddo ! i and k loops
if (associated(ADp%du_dt_str)) then
do i=is,ie ; if (abs(ADp%du_dt_str(I,j,nz)) < accel_underflow) ADp%du_dt_str(I,j,nz) = 0.0 ; enddo
do k=nz-1,1,-1 ; do I=Isq,Ieq ; if (do_i(I)) then
ADp%du_dt_str(I,j,k) = ADp%du_dt_str(I,j,k) + c1(I,k+1) * ADp%du_dt_str(I,j,k+1)
if (abs(ADp%du_dt_str(I,j,k)) < accel_underflow) ADp%du_dt_str(I,j,k) = 0.0
endif ; enddo ; enddo
endif
if (associated(ADp%du_dt_visc)) then ; do k=1,nz ; do I=Isq,Ieq
ADp%du_dt_visc(I,j,k) = (u(I,j,k) - ADp%du_dt_visc(I,j,k))*Idt
if (abs(ADp%du_dt_visc(I,j,k)) < accel_underflow) ADp%du_dt_visc(I,j,k) = 0.0
enddo ; enddo ; endif
if (associated(visc%taux_shelf)) then ; do I=Isq,Ieq
visc%taux_shelf(I,j) = -GV%Rho0*CS%a1_shelf_u(I,j)*u(I,j,1) ! - u_shelf?
enddo ; endif
if (PRESENT(taux_bot)) then
do I=Isq,Ieq
taux_bot(I,j) = GV%Rho0 * (u(I,j,nz)*CS%a_u(I,j,nz+1))
enddo
if (CS%Channel_drag) then ; do k=1,nz ; do I=Isq,Ieq
taux_bot(I,j) = taux_bot(I,j) + GV%Rho0 * (Ray(I,k)*u(I,j,k))
enddo ; enddo ; endif
endif
! When mixing down Eulerian current + Stokes drift subtract after calling solver
if (DoStokesMixing) then ; do k=1,nz ; do I=Isq,Ieq
if (do_i(I)) u(I,j,k) = u(I,j,k) - Waves%Us_x(I,j,k)
enddo ; enddo ; endif
enddo ! end u-component j loop
! Now work on the meridional velocity component.
!$OMP parallel do default(shared) firstprivate(Ray) &
!$OMP private(do_i,surface_stress,zDS,stress,h_a,hfr, &
!$OMP b_denom_1,b1,d1,c1)
do J=Jsq,Jeq
do i=is,ie ; do_i(i) = (G%mask2dCv(i,J) > 0) ; enddo
! When mixing down Eulerian current + Stokes drift add before calling solver
if (DoStokesMixing) then ; do k=1,nz ; do i=is,ie
if (do_i(i)) v(i,j,k) = v(i,j,k) + Waves%Us_y(i,j,k)
enddo ; enddo ; endif
if (associated(ADp%dv_dt_visc)) then ; do k=1,nz ; do i=is,ie
ADp%dv_dt_visc(i,J,k) = v(i,J,k)
enddo ; enddo ; endif
if (associated(ADp%dv_dt_str)) then ; do k=1,nz ; do i=is,ie
ADp%dv_dt_str(i,J,k) = 0.0
enddo ; enddo ; endif
! One option is to have the wind stress applied as a body force
! over the topmost Hmix fluid. If DIRECT_STRESS is not defined,
! the wind stress is applied as a stress boundary condition.
if (CS%direct_stress) then
do i=is,ie ; if (do_i(i)) then
surface_stress(i) = 0.0
zDS = 0.0
stress = dt_Rho0 * forces%tauy(i,J)
do k=1,nz
h_a = 0.5 * (h(i,J,k) + h(i,J+1,k)) + h_neglect
hfr = 1.0 ; if ((zDS+h_a) > Hmix) hfr = (Hmix - zDS) / h_a
v(i,J,k) = v(i,J,k) + I_Hmix * hfr * stress
if (associated(ADp%dv_dt_str)) ADp%dv_dt_str(i,J,k) = (I_Hmix * hfr * stress) * Idt
zDS = zDS + h_a ; if (zDS >= Hmix) exit
enddo
endif ; enddo ! end of i loop
else ; do i=is,ie
surface_stress(i) = dt_Rho0 * (G%mask2dCv(i,J)*forces%tauy(i,J))
enddo ; endif ! direct_stress
if (CS%Channel_drag) then ; do k=1,nz ; do i=is,ie
Ray(i,k) = visc%Ray_v(i,J,k)
enddo ; enddo ; endif
do i=is,ie ; if (do_i(i)) then
b_denom_1 = CS%h_v(i,J,1) + dt_Z_to_H * (Ray(i,1) + CS%a_v(i,J,1))
b1(i) = 1.0 / (b_denom_1 + dt_Z_to_H*CS%a_v(i,J,2))
d1(i) = b_denom_1 * b1(i)
v(i,J,1) = b1(i) * (CS%h_v(i,J,1) * v(i,J,1) + surface_stress(i))
if (associated(ADp%dv_dt_str)) &
ADp%dv_dt_str(i,J,1) = b1(i) * (CS%h_v(i,J,1) * ADp%dv_dt_str(i,J,1) + surface_stress(i)*Idt)
endif ; enddo
do k=2,nz ; do i=is,ie ; if (do_i(i)) then
c1(i,k) = dt_Z_to_H * CS%a_v(i,J,K) * b1(i)
b_denom_1 = CS%h_v(i,J,k) + dt_Z_to_H * (Ray(i,k) + CS%a_v(i,J,K)*d1(i))
b1(i) = 1.0 / (b_denom_1 + dt_Z_to_H * CS%a_v(i,J,K+1))
d1(i) = b_denom_1 * b1(i)
v(i,J,k) = (CS%h_v(i,J,k) * v(i,J,k) + dt_Z_to_H * CS%a_v(i,J,K) * v(i,J,k-1)) * b1(i)
if (associated(ADp%dv_dt_str)) &
ADp%dv_dt_str(i,J,k) = (CS%h_v(i,J,k) * ADp%dv_dt_str(i,J,k) + &
dt_Z_to_H * CS%a_v(i,J,K) * ADp%dv_dt_str(i,J,k-1)) * b1(i)
endif ; enddo ; enddo
do k=nz-1,1,-1 ; do i=is,ie ; if (do_i(i)) then
v(i,J,k) = v(i,J,k) + c1(i,k+1) * v(i,J,k+1)
endif ; enddo ; enddo ! i and k loops
if (associated(ADp%dv_dt_str)) then
do i=is,ie ; if (abs(ADp%dv_dt_str(i,J,nz)) < accel_underflow) ADp%dv_dt_str(i,J,nz) = 0.0 ; enddo
do k=nz-1,1,-1 ; do i=is,ie ; if (do_i(i)) then
ADp%dv_dt_str(i,J,k) = ADp%dv_dt_str(i,J,k) + c1(i,k+1) * ADp%dv_dt_str(i,J,k+1)
if (abs(ADp%dv_dt_str(i,J,k)) < accel_underflow) ADp%dv_dt_str(i,J,k) = 0.0
endif ; enddo ; enddo
endif
if (associated(ADp%dv_dt_visc)) then ; do k=1,nz ; do i=is,ie
ADp%dv_dt_visc(i,J,k) = (v(i,J,k) - ADp%dv_dt_visc(i,J,k))*Idt
if (abs(ADp%dv_dt_visc(i,J,k)) < accel_underflow) ADp%dv_dt_visc(i,J,k) = 0.0
enddo ; enddo ; endif
if (associated(visc%tauy_shelf)) then ; do i=is,ie
visc%tauy_shelf(i,J) = -GV%Rho0*CS%a1_shelf_v(i,J)*v(i,J,1) ! - v_shelf?
enddo ; endif
if (present(tauy_bot)) then
do i=is,ie
tauy_bot(i,J) = GV%Rho0 * (v(i,J,nz)*CS%a_v(i,J,nz+1))
enddo
if (CS%Channel_drag) then ; do k=1,nz ; do i=is,ie
tauy_bot(i,J) = tauy_bot(i,J) + GV%Rho0 * (Ray(i,k)*v(i,J,k))
enddo ; enddo ; endif
endif
! When mixing down Eulerian current + Stokes drift subtract after calling solver
if (DoStokesMixing) then ; do k=1,nz ; do i=is,ie
if (do_i(i)) v(i,J,k) = v(i,J,k) - Waves%Us_y(i,J,k)
enddo ; enddo ; endif
enddo ! end of v-component J loop
call vertvisc_limit_vel(u, v, h, ADp, CDp, forces, visc, dt, G, GV, US, CS)
! Here the velocities associated with open boundary conditions are applied.
if (associated(OBC)) then
do n=1,OBC%number_of_segments
if (OBC%segment(n)%specified) then
if (OBC%segment(n)%is_N_or_S) then
J = OBC%segment(n)%HI%JsdB
do k=1,nz ; do i=OBC%segment(n)%HI%isd,OBC%segment(n)%HI%ied
v(i,J,k) = OBC%segment(n)%normal_vel(i,J,k)
enddo ; enddo
elseif (OBC%segment(n)%is_E_or_W) then
I = OBC%segment(n)%HI%IsdB
do k=1,nz ; do j=OBC%segment(n)%HI%jsd,OBC%segment(n)%HI%jed
u(I,j,k) = OBC%segment(n)%normal_vel(I,j,k)
enddo ; enddo
endif
endif
enddo
endif
! Offer diagnostic fields for averaging.
if (CS%id_du_dt_visc > 0) &
call post_data(CS%id_du_dt_visc, ADp%du_dt_visc, CS%diag)
if (CS%id_dv_dt_visc > 0) &
call post_data(CS%id_dv_dt_visc, ADp%dv_dt_visc, CS%diag)
if (present(taux_bot) .and. (CS%id_taux_bot > 0)) &
call post_data(CS%id_taux_bot, taux_bot, CS%diag)
if (present(tauy_bot) .and. (CS%id_tauy_bot > 0)) &
call post_data(CS%id_tauy_bot, tauy_bot, CS%diag)
if (CS%id_du_dt_str > 0) &
call post_data(CS%id_du_dt_str, ADp%du_dt_str, CS%diag)
if (CS%id_dv_dt_str > 0) &
call post_data(CS%id_dv_dt_str, ADp%dv_dt_str, CS%diag)
if (associated(ADp%du_dt_visc) .and. associated(ADp%du_dt_visc)) then
! Diagnostics of the fractional thicknesses times momentum budget terms
! 3D diagnostics of hf_du(dv)_dt_visc are commented because there is no clarity on proper remapping grid option.
! The code is retained for debugging purposes in the future.
!if (CS%id_hf_du_dt_visc > 0) &
! call post_product_u(CS%id_hf_du_dt_visc, ADp%du_dt_visc, ADp%diag_hfrac_u, G, nz, CS%diag)
!if (CS%id_hf_dv_dt_visc > 0) &
! call post_product_v(CS%id_hf_dv_dt_visc, ADp%dv_dt_visc, ADp%diag_hfrac_v, G, nz, CS%diag)
! Diagnostics for thickness-weighted vertically averaged viscous accelerations
if (CS%id_hf_du_dt_visc_2d > 0) &
call post_product_sum_u(CS%id_hf_du_dt_visc_2d, ADp%du_dt_visc, ADp%diag_hfrac_u, G, nz, CS%diag)
if (CS%id_hf_dv_dt_visc_2d > 0) &
call post_product_sum_v(CS%id_hf_dv_dt_visc_2d, ADp%dv_dt_visc, ADp%diag_hfrac_v, G, nz, CS%diag)
! Diagnostics for thickness x viscous accelerations
if (CS%id_h_du_dt_visc > 0) call post_product_u(CS%id_h_du_dt_visc, ADp%du_dt_visc, ADp%diag_hu, G, nz, CS%diag)
if (CS%id_h_dv_dt_visc > 0) call post_product_v(CS%id_h_dv_dt_visc, ADp%dv_dt_visc, ADp%diag_hv, G, nz, CS%diag)
endif
if (associated(ADp%du_dt_str) .and. associated(ADp%dv_dt_str)) then
! Diagnostics for thickness x wind stress accelerations
if (CS%id_h_du_dt_str > 0) call post_product_u(CS%id_h_du_dt_str, ADp%du_dt_str, ADp%diag_hu, G, nz, CS%diag)
if (CS%id_h_dv_dt_str > 0) call post_product_v(CS%id_h_dv_dt_str, ADp%dv_dt_str, ADp%diag_hv, G, nz, CS%diag)
! Diagnostics for wind stress accelerations multiplied by visc_rem_[uv],
if (CS%id_du_dt_str_visc_rem > 0) &
call post_product_u(CS%id_du_dt_str_visc_rem, ADp%du_dt_str, ADp%visc_rem_u, G, nz, CS%diag)
if (CS%id_dv_dt_str_visc_rem > 0) &
call post_product_v(CS%id_dv_dt_str_visc_rem, ADp%dv_dt_str, ADp%visc_rem_v, G, nz, CS%diag)
endif
end subroutine vertvisc
!> Calculate the fraction of momentum originally in a layer that remains in the water column
!! after a time-step of viscosity, equivalently the fraction of a time-step's worth of
!! barotropic acceleration that a layer experiences after viscosity is applied.
subroutine vertvisc_remnant(visc, visc_rem_u, visc_rem_v, dt, G, GV, US, CS)
type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
type(vertvisc_type), intent(in) :: visc !< Viscosities and bottom drag
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
intent(inout) :: visc_rem_u !< Fraction of a time-step's worth of a
!! barotropic acceleration that a layer experiences after
!! viscosity is applied in the zonal direction [nondim]
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
intent(inout) :: visc_rem_v !< Fraction of a time-step's worth of a
!! barotropic acceleration that a layer experiences after
!! viscosity is applied in the meridional direction [nondim]
real, intent(in) :: dt !< Time increment [T ~> s]
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
type(vertvisc_CS), pointer :: CS !< Vertical viscosity control structure
! Local variables
real :: b1(SZIB_(G)) ! A variable used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1].
real :: c1(SZIB_(G),SZK_(GV)) ! A variable used by the tridiagonal solver [nondim].
real :: d1(SZIB_(G)) ! d1=1-c1 is used by the tridiagonal solver [nondim].
real :: Ray(SZIB_(G),SZK_(GV)) ! Ray is the Rayleigh-drag velocity [Z T-1 ~> m s-1].
real :: b_denom_1 ! The first term in the denominator of b1 [H ~> m or kg m-2].
real :: dt_Z_to_H ! The time step times the conversion from Z to the
! units of thickness [T H Z-1 ~> s or s kg m-3].
logical :: do_i(SZIB_(G))
integer :: i, j, k, is, ie, Isq, Ieq, Jsq, Jeq, nz
is = G%isc ; ie = G%iec
Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB ; nz = GV%ke
if (.not.associated(CS)) call MOM_error(FATAL,"MOM_vert_friction(visc): "// &
"Module must be initialized before it is used.")
if (.not.CS%initialized) call MOM_error(FATAL,"MOM_vert_friction(remant): "// &
"Module must be initialized before it is used.")
dt_Z_to_H = dt*GV%Z_to_H
do k=1,nz ; do i=Isq,Ieq ; Ray(i,k) = 0.0 ; enddo ; enddo
! Find the zonal viscous remnant using a modification of a standard tridagonal solver.
!$OMP parallel do default(shared) firstprivate(Ray) private(do_i,b_denom_1,b1,d1,c1)
do j=G%jsc,G%jec
do I=Isq,Ieq ; do_i(I) = (G%mask2dCu(I,j) > 0) ; enddo
if (CS%Channel_drag) then ; do k=1,nz ; do I=Isq,Ieq
Ray(I,k) = visc%Ray_u(I,j,k)
enddo ; enddo ; endif
do I=Isq,Ieq ; if (do_i(I)) then
b_denom_1 = CS%h_u(I,j,1) + dt_Z_to_H * (Ray(I,1) + CS%a_u(I,j,1))
b1(I) = 1.0 / (b_denom_1 + dt_Z_to_H*CS%a_u(I,j,2))
d1(I) = b_denom_1 * b1(I)
visc_rem_u(I,j,1) = b1(I) * CS%h_u(I,j,1)
endif ; enddo
do k=2,nz ; do I=Isq,Ieq ; if (do_i(I)) then
c1(I,k) = dt_Z_to_H * CS%a_u(I,j,K)*b1(I)
b_denom_1 = CS%h_u(I,j,k) + dt_Z_to_H * (Ray(I,k) + CS%a_u(I,j,K)*d1(I))
b1(I) = 1.0 / (b_denom_1 + dt_Z_to_H * CS%a_u(I,j,K+1))
d1(I) = b_denom_1 * b1(I)
visc_rem_u(I,j,k) = (CS%h_u(I,j,k) + dt_Z_to_H * CS%a_u(I,j,K) * visc_rem_u(I,j,k-1)) * b1(I)
endif ; enddo ; enddo
do k=nz-1,1,-1 ; do I=Isq,Ieq ; if (do_i(I)) then
visc_rem_u(I,j,k) = visc_rem_u(I,j,k) + c1(I,k+1)*visc_rem_u(I,j,k+1)
endif ; enddo ; enddo ! i and k loops
enddo ! end u-component j loop
! Now find the meridional viscous remnant using the robust tridiagonal solver.
!$OMP parallel do default(shared) firstprivate(Ray) private(do_i,b_denom_1,b1,d1,c1)
do J=Jsq,Jeq
do i=is,ie ; do_i(i) = (G%mask2dCv(i,J) > 0) ; enddo
if (CS%Channel_drag) then ; do k=1,nz ; do i=is,ie
Ray(i,k) = visc%Ray_v(i,J,k)
enddo ; enddo ; endif
do i=is,ie ; if (do_i(i)) then
b_denom_1 = CS%h_v(i,J,1) + dt_Z_to_H * (Ray(i,1) + CS%a_v(i,J,1))
b1(i) = 1.0 / (b_denom_1 + dt_Z_to_H*CS%a_v(i,J,2))
d1(i) = b_denom_1 * b1(i)
visc_rem_v(i,J,1) = b1(i) * CS%h_v(i,J,1)
endif ; enddo
do k=2,nz ; do i=is,ie ; if (do_i(i)) then
c1(i,k) = dt_Z_to_H * CS%a_v(i,J,K)*b1(i)
b_denom_1 = CS%h_v(i,J,k) + dt_Z_to_H * (Ray(i,k) + CS%a_v(i,J,K)*d1(i))
b1(i) = 1.0 / (b_denom_1 + dt_Z_to_H * CS%a_v(i,J,K+1))
d1(i) = b_denom_1 * b1(i)
visc_rem_v(i,J,k) = (CS%h_v(i,J,k) + dt_Z_to_H * CS%a_v(i,J,K) * visc_rem_v(i,J,k-1)) * b1(i)
endif ; enddo ; enddo
do k=nz-1,1,-1 ; do i=is,ie ; if (do_i(i)) then
visc_rem_v(i,J,k) = visc_rem_v(i,J,k) + c1(i,k+1)*visc_rem_v(i,J,k+1)
endif ; enddo ; enddo ! i and k loops
enddo ! end of v-component J loop
if (CS%debug) then
call uvchksum("visc_rem_[uv]", visc_rem_u, visc_rem_v, G%HI, haloshift=0, &
scalar_pair=.true.)
endif
end subroutine vertvisc_remnant
!> Calculate the coupling coefficients (CS%a_u and CS%a_v)
!! and effective layer thicknesses (CS%h_u and CS%h_v) for later use in the
!! applying the implicit vertical viscosity via vertvisc().
subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC)
type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1]
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
type(vertvisc_type), intent(in) :: visc !< Viscosities and bottom drag
real, intent(in) :: dt !< Time increment [T ~> s]
type(vertvisc_CS), pointer :: CS !< Vertical viscosity control structure
type(ocean_OBC_type), pointer :: OBC !< Open boundary condition structure
! Field from forces used in this subroutine:
! ustar: the friction velocity [Z T-1 ~> m s-1], used here as the mixing
! velocity in the mixed layer if NKML > 1 in a bulk mixed layer.
! Local variables
real, dimension(SZIB_(G),SZK_(GV)) :: &
h_harm, & ! Harmonic mean of the thicknesses around a velocity grid point,
! given by 2*(h+ * h-)/(h+ + h-) [H ~> m or kg m-2].
h_arith, & ! The arithmetic mean thickness [H ~> m or kg m-2].
h_delta, & ! The lateral difference of thickness [H ~> m or kg m-2].
hvel, & ! hvel is the thickness used at a velocity grid point [H ~> m or kg m-2].
hvel_shelf ! The equivalent of hvel under shelves [H ~> m or kg m-2].
real, dimension(SZIB_(G),SZK_(GV)+1) :: &
a_cpl, & ! The drag coefficients across interfaces [Z T-1 ~> m s-1]. a_cpl times
! the velocity difference gives the stress across an interface.
a_shelf, & ! The drag coefficients across interfaces in water columns under
! ice shelves [Z T-1 ~> m s-1].
z_i ! An estimate of each interface's height above the bottom,
! normalized by the bottom boundary layer thickness [nondim]
real, dimension(SZIB_(G)) :: &
kv_bbl, & ! The bottom boundary layer viscosity [Z2 T-1 ~> m2 s-1].
bbl_thick, & ! The bottom boundary layer thickness [H ~> m or kg m-2].
I_Hbbl, & ! The inverse of the bottom boundary layer thickness [H-1 ~> m-1 or m2 kg-1].
I_Htbl, & ! The inverse of the top boundary layer thickness [H-1 ~> m-1 or m2 kg-1].
zcol1, & ! The height of the interfaces to the north and south of a
zcol2, & ! v-point [H ~> m or kg m-2].
Ztop_min, & ! The deeper of the two adjacent surface heights [H ~> m or kg m-2].
Dmin, & ! The shallower of the two adjacent bottom depths converted to
! thickness units [H ~> m or kg m-2].
zh, & ! An estimate of the interface's distance from the bottom
! based on harmonic mean thicknesses [H ~> m or kg m-2].
h_ml ! The mixed layer depth [H ~> m or kg m-2].
real, allocatable, dimension(:,:) :: hML_u ! Diagnostic of the mixed layer depth at u points [H ~> m or kg m-2].
real, allocatable, dimension(:,:) :: hML_v ! Diagnostic of the mixed layer depth at v points [H ~> m or kg m-2].
real, allocatable, dimension(:,:,:) :: Kv_u !< Total vertical viscosity at u-points [Z2 T-1 ~> m2 s-1].
real, allocatable, dimension(:,:,:) :: Kv_v !< Total vertical viscosity at v-points [Z2 T-1 ~> m2 s-1].
real :: zcol(SZI_(G)) ! The height of an interface at h-points [H ~> m or kg m-2].
real :: botfn ! A function which goes from 1 at the bottom to 0 much more
! than Hbbl into the interior.
real :: topfn ! A function which goes from 1 at the top to 0 much more
! than Htbl into the interior.
real :: z2 ! The distance from the bottom, normalized by Hbbl [nondim]
real :: z2_wt ! A nondimensional (0-1) weight used when calculating z2.
real :: z_clear ! The clearance of an interface above the surrounding topography [H ~> m or kg m-2].
real :: a_cpl_max ! The maximum drag coefficient across interfaces, set so that it will be
! representable as a 32-bit float in MKS units [Z T-1 ~> m s-1]
real :: h_neglect ! A thickness that is so small it is usually lost
! in roundoff and can be neglected [H ~> m or kg m-2].
real :: I_valBL ! The inverse of a scaling factor determining when water is
! still within the boundary layer, as determined by the sum
! of the harmonic mean thicknesses.
logical, dimension(SZIB_(G)) :: do_i, do_i_shelf
logical :: do_any_shelf
integer, dimension(SZIB_(G)) :: &
zi_dir ! A trinary logical array indicating which thicknesses to use for
! finding z_clear.
integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB ; nz = GV%ke
if (.not.associated(CS)) call MOM_error(FATAL,"MOM_vert_friction(coef): "// &
"Module must be initialized before it is used.")
if (.not.CS%initialized) call MOM_error(FATAL,"MOM_vert_friction(coef): "// &
"Module must be initialized before it is used.")
h_neglect = GV%H_subroundoff
a_cpl_max = 1.0e37 * US%m_to_Z * US%T_to_s
I_Hbbl(:) = 1.0 / (CS%Hbbl + h_neglect)
I_valBL = 0.0 ; if (CS%harm_BL_val > 0.0) I_valBL = 1.0 / CS%harm_BL_val
if (CS%id_Kv_u > 0) allocate(Kv_u(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke), source=0.0)
if (CS%id_Kv_v > 0) allocate(Kv_v(G%isd:G%ied,G%JsdB:G%JedB,GV%ke), source=0.0)
if (CS%debug .or. (CS%id_hML_u > 0)) allocate(hML_u(G%IsdB:G%IedB,G%jsd:G%jed), source=0.0)
if (CS%debug .or. (CS%id_hML_v > 0)) allocate(hML_v(G%isd:G%ied,G%JsdB:G%JedB), source=0.0)
if ((associated(visc%taux_shelf) .or. associated(forces%frac_shelf_u)) .and. &
.not.associated(CS%a1_shelf_u)) then
allocate(CS%a1_shelf_u(G%IsdB:G%IedB,G%jsd:G%jed), source=0.0)
endif
if ((associated(visc%tauy_shelf) .or. associated(forces%frac_shelf_v)) .and. &
.not.associated(CS%a1_shelf_v)) then
allocate(CS%a1_shelf_v(G%isd:G%ied,G%JsdB:G%JedB), source=0.0)
endif
!$OMP parallel do default(private) shared(G,GV,US,CS,visc,Isq,Ieq,nz,u,h,forces,hML_u, &
!$OMP OBC,h_neglect,dt,I_valBL,Kv_u,a_cpl_max) &
!$OMP firstprivate(i_hbbl)
do j=G%Jsc,G%Jec
do I=Isq,Ieq ; do_i(I) = (G%mask2dCu(I,j) > 0) ; enddo
if (CS%bottomdraglaw) then ; do I=Isq,Ieq
kv_bbl(I) = visc%Kv_bbl_u(I,j)
bbl_thick(I) = visc%bbl_thick_u(I,j) * GV%Z_to_H + h_neglect
if (do_i(I)) I_Hbbl(I) = 1.0 / bbl_thick(I)
enddo ; endif
do k=1,nz ; do I=Isq,Ieq ; if (do_i(I)) then
h_harm(I,k) = 2.0*h(i,j,k)*h(i+1,j,k) / (h(i,j,k)+h(i+1,j,k)+h_neglect)
h_arith(I,k) = 0.5*(h(i+1,j,k)+h(i,j,k))
h_delta(I,k) = h(i+1,j,k) - h(i,j,k)
endif ; enddo ; enddo
do I=Isq,Ieq
Dmin(I) = min(G%bathyT(i,j), G%bathyT(i+1,j)) * GV%Z_to_H
zi_dir(I) = 0
enddo
! Project thickness outward across OBCs using a zero-gradient condition.
if (associated(OBC)) then ; if (OBC%number_of_segments > 0) then
do I=Isq,Ieq ; if (do_i(I) .and. (OBC%segnum_u(I,j) /= OBC_NONE)) then
if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_E) then
do k=1,nz ; h_harm(I,k) = h(i,j,k) ; h_arith(I,k) = h(i,j,k) ; h_delta(I,k) = 0. ; enddo
Dmin(I) = G%bathyT(i,j) * GV%Z_to_H
zi_dir(I) = -1
elseif (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_W) then
do k=1,nz ; h_harm(I,k) = h(i+1,j,k) ; h_arith(I,k) = h(i+1,j,k) ; h_delta(I,k) = 0. ; enddo
Dmin(I) = G%bathyT(i+1,j) * GV%Z_to_H
zi_dir(I) = 1
endif
endif ; enddo
endif ; endif
! The following block calculates the thicknesses at velocity
! grid points for the vertical viscosity (hvel). Near the
! bottom an upwind biased thickness is used to control the effect
! of spurious Montgomery potential gradients at the bottom where
! nearly massless layers layers ride over the topography.
if (CS%harmonic_visc) then
do I=Isq,Ieq ; z_i(I,nz+1) = 0.0 ; enddo
do k=nz,1,-1 ; do I=Isq,Ieq ; if (do_i(I)) then
hvel(I,k) = h_harm(I,k)
if (u(I,j,k) * h_delta(I,k) < 0) then
z2 = z_i(I,k+1) ; botfn = 1.0 / (1.0 + 0.09*z2*z2*z2*z2*z2*z2)
hvel(I,k) = (1.0-botfn)*h_harm(I,k) + botfn*h_arith(I,k)
endif
z_i(I,k) = z_i(I,k+1) + h_harm(I,k)*I_Hbbl(I)
endif ; enddo ; enddo ! i & k loops
else ! Not harmonic_visc
do I=Isq,Ieq ; zh(I) = 0.0 ; z_i(I,nz+1) = 0.0 ; enddo
do i=Isq,Ieq+1 ; zcol(i) = -G%bathyT(i,j) * GV%Z_to_H ; enddo
do k=nz,1,-1
do i=Isq,Ieq+1 ; zcol(i) = zcol(i) + h(i,j,k) ; enddo
do I=Isq,Ieq ; if (do_i(I)) then
zh(I) = zh(I) + h_harm(I,k)
z_clear = max(zcol(i),zcol(i+1)) + Dmin(I)
if (zi_dir(I) < 0) z_clear = zcol(i) + Dmin(I)
if (zi_dir(I) > 0) z_clear = zcol(i+1) + Dmin(I)
z_i(I,k) = max(zh(I), z_clear) * I_Hbbl(I)
hvel(I,k) = h_arith(I,k)
if (u(I,j,k) * h_delta(I,k) > 0) then
if (zh(I) * I_Hbbl(I) < CS%harm_BL_val) then
hvel(I,k) = h_harm(I,k)
else
z2_wt = 1.0 ; if (zh(I) * I_Hbbl(I) < 2.0*CS%harm_BL_val) &
z2_wt = max(0.0, min(1.0, zh(I) * I_Hbbl(I) * I_valBL - 1.0))
z2 = z2_wt * (max(zh(I), z_clear) * I_Hbbl(I))
botfn = 1.0 / (1.0 + 0.09*z2*z2*z2*z2*z2*z2)
hvel(I,k) = (1.0-botfn)*h_arith(I,k) + botfn*h_harm(I,k)
endif
endif
endif ; enddo ! i loop
enddo ! k loop
endif
call find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, h_ml, &
dt, j, G, GV, US, CS, visc, forces, work_on_u=.true., OBC=OBC)
if (allocated(hML_u)) then
do i=isq,ieq ; if (do_i(i)) then ; hML_u(I,j) = h_ml(I) ; endif ; enddo
endif
do_any_shelf = .false.
if (associated(forces%frac_shelf_u)) then
do I=Isq,Ieq
CS%a1_shelf_u(I,j) = 0.0
do_i_shelf(I) = (do_i(I) .and. forces%frac_shelf_u(I,j) > 0.0)
if (do_i_shelf(I)) do_any_shelf = .true.
enddo
if (do_any_shelf) then
if (CS%harmonic_visc) then
do k=1,nz ; do I=Isq,Ieq ; hvel_shelf(I,k) = hvel(I,k) ; enddo ; enddo
else ! Find upwind-biased thickness near the surface.
! Perhaps this needs to be done more carefully, via find_eta.
do I=Isq,Ieq ; if (do_i_shelf(I)) then
zh(I) = 0.0 ; Ztop_min(I) = min(zcol(i), zcol(i+1))
I_HTbl(I) = 1.0 / (visc%tbl_thick_shelf_u(I,j)*GV%Z_to_H + h_neglect)
endif ; enddo
do k=1,nz
do i=Isq,Ieq+1 ; zcol(i) = zcol(i) - h(i,j,k) ; enddo
do I=Isq,Ieq ; if (do_i_shelf(I)) then
zh(I) = zh(I) + h_harm(I,k)
hvel_shelf(I,k) = hvel(I,k)
if (u(I,j,k) * h_delta(I,k) > 0) then
if (zh(I) * I_HTbl(I) < CS%harm_BL_val) then
hvel_shelf(I,k) = min(hvel(I,k), h_harm(I,k))
else
z2_wt = 1.0 ; if (zh(I) * I_HTbl(I) < 2.0*CS%harm_BL_val) &
z2_wt = max(0.0, min(1.0, zh(I) * I_HTbl(I) * I_valBL - 1.0))
z2 = z2_wt * (max(zh(I), Ztop_min(I) - min(zcol(i),zcol(i+1))) * I_HTbl(I))
topfn = 1.0 / (1.0 + 0.09*z2**6)
hvel_shelf(I,k) = min(hvel(I,k), (1.0-topfn)*h_arith(I,k) + topfn*h_harm(I,k))
endif
endif
endif ; enddo
enddo
endif
call find_coupling_coef(a_shelf, hvel_shelf, do_i_shelf, h_harm, bbl_thick, &
kv_bbl, z_i, h_ml, dt, j, G, GV, US, CS, visc, forces, &
work_on_u=.true., OBC=OBC, shelf=.true.)
do I=Isq,Ieq ; if (do_i_shelf(I)) CS%a1_shelf_u(I,j) = a_shelf(I,1) ; enddo
endif
endif
if (do_any_shelf) then
do K=1,nz+1 ; do I=Isq,Ieq ; if (do_i_shelf(I)) then
CS%a_u(I,j,K) = min(a_cpl_max, forces%frac_shelf_u(I,j) * a_shelf(I,K) + &
(1.0-forces%frac_shelf_u(I,j)) * a_cpl(I,K))
! This is Alistair's suggestion, but it destabilizes the model. I do not know why. RWH
! CS%a_u(I,j,K) = min(a_cpl_max, forces%frac_shelf_u(I,j) * max(a_shelf(I,K), a_cpl(I,K)) + &
! (1.0-forces%frac_shelf_u(I,j)) * a_cpl(I,K))
elseif (do_i(I)) then
CS%a_u(I,j,K) = min(a_cpl_max, a_cpl(I,K))
endif ; enddo ; enddo
do k=1,nz ; do I=Isq,Ieq ; if (do_i_shelf(I)) then
! Should we instead take the inverse of the average of the inverses?
CS%h_u(I,j,k) = forces%frac_shelf_u(I,j) * hvel_shelf(I,k) + &
(1.0-forces%frac_shelf_u(I,j)) * hvel(I,k) + h_neglect
elseif (do_i(I)) then
CS%h_u(I,j,k) = hvel(I,k) + h_neglect
endif ; enddo ; enddo
else
do K=1,nz+1 ; do I=Isq,Ieq ; if (do_i(I)) CS%a_u(I,j,K) = min(a_cpl_max, a_cpl(I,K)) ; enddo ; enddo
do k=1,nz ; do I=Isq,Ieq ; if (do_i(I)) CS%h_u(I,j,k) = hvel(I,k) + h_neglect ; enddo ; enddo
endif
! Diagnose total Kv at u-points
if (CS%id_Kv_u > 0) then
do k=1,nz ; do I=Isq,Ieq
if (do_i(I)) Kv_u(I,j,k) = 0.5 * GV%H_to_Z*(CS%a_u(I,j,K)+CS%a_u(I,j,K+1)) * CS%h_u(I,j,k)
enddo ; enddo
endif
enddo
! Now work on v-points.
!$OMP parallel do default(private) shared(G,GV,CS,US,visc,is,ie,Jsq,Jeq,nz,v,h,forces,hML_v, &
!$OMP OBC,h_neglect,dt,I_valBL,Kv_v,a_cpl_max) &
!$OMP firstprivate(i_hbbl)
do J=Jsq,Jeq
do i=is,ie ; do_i(i) = (G%mask2dCv(i,J) > 0) ; enddo
if (CS%bottomdraglaw) then ; do i=is,ie
kv_bbl(i) = visc%Kv_bbl_v(i,J)
bbl_thick(i) = visc%bbl_thick_v(i,J) * GV%Z_to_H + h_neglect
if (do_i(i)) I_Hbbl(i) = 1.0 / bbl_thick(i)
enddo ; endif
do k=1,nz ; do i=is,ie ; if (do_i(i)) then
h_harm(i,k) = 2.0*h(i,j,k)*h(i,j+1,k) / (h(i,j,k)+h(i,j+1,k)+h_neglect)
h_arith(i,k) = 0.5*(h(i,j+1,k)+h(i,j,k))
h_delta(i,k) = h(i,j+1,k) - h(i,j,k)
endif ; enddo ; enddo
do i=is,ie
Dmin(i) = min(G%bathyT(i,j), G%bathyT(i,j+1)) * GV%Z_to_H
zi_dir(i) = 0
enddo
! Project thickness outward across OBCs using a zero-gradient condition.
if (associated(OBC)) then ; if (OBC%number_of_segments > 0) then
do i=is,ie ; if (do_i(i) .and. (OBC%segnum_v(i,J) /= OBC_NONE)) then
if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_N) then
do k=1,nz ; h_harm(I,k) = h(i,j,k) ; h_arith(I,k) = h(i,j,k) ; h_delta(i,k) = 0. ; enddo
Dmin(I) = G%bathyT(i,j) * GV%Z_to_H
zi_dir(I) = -1
elseif (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_S) then
do k=1,nz ; h_harm(i,k) = h(i,j+1,k) ; h_arith(i,k) = h(i,j+1,k) ; h_delta(i,k) = 0. ; enddo
Dmin(i) = G%bathyT(i,j+1) * GV%Z_to_H
zi_dir(i) = 1
endif
endif ; enddo
endif ; endif
! The following block calculates the thicknesses at velocity
! grid points for the vertical viscosity (hvel). Near the
! bottom an upwind biased thickness is used to control the effect
! of spurious Montgomery potential gradients at the bottom where
! nearly massless layers layers ride over the topography.
if (CS%harmonic_visc) then
do i=is,ie ; z_i(i,nz+1) = 0.0 ; enddo
do k=nz,1,-1 ; do i=is,ie ; if (do_i(i)) then
hvel(i,k) = h_harm(i,k)
if (v(i,J,k) * h_delta(i,k) < 0) then
z2 = z_i(i,k+1) ; botfn = 1.0 / (1.0 + 0.09*z2*z2*z2*z2*z2*z2)
hvel(i,k) = (1.0-botfn)*h_harm(i,k) + botfn*h_arith(i,k)
endif
z_i(i,k) = z_i(i,k+1) + h_harm(i,k)*I_Hbbl(i)
endif ; enddo ; enddo ! i & k loops
else ! Not harmonic_visc
do i=is,ie
zh(i) = 0.0 ; z_i(i,nz+1) = 0.0
zcol1(i) = -G%bathyT(i,j) * GV%Z_to_H
zcol2(i) = -G%bathyT(i,j+1) * GV%Z_to_H
enddo
do k=nz,1,-1 ; do i=is,ie ; if (do_i(i)) then
zh(i) = zh(i) + h_harm(i,k)
zcol1(i) = zcol1(i) + h(i,j,k) ; zcol2(i) = zcol2(i) + h(i,j+1,k)
z_clear = max(zcol1(i),zcol2(i)) + Dmin(i)
if (zi_dir(i) < 0) z_clear = zcol1(i) + Dmin(I)
if (zi_dir(i) > 0) z_clear = zcol2(i) + Dmin(I)
z_i(I,k) = max(zh(i), z_clear) * I_Hbbl(i)
hvel(i,k) = h_arith(i,k)
if (v(i,J,k) * h_delta(i,k) > 0) then