forked from mom-ocean/MOM6
-
Notifications
You must be signed in to change notification settings - Fork 51
/
MOM_kappa_shear.F90
2078 lines (1893 loc) · 107 KB
/
MOM_kappa_shear.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
!> Shear-dependent mixing following Jackson et al. 2008.
module MOM_kappa_shear
! This file is part of MOM6. See LICENSE.md for the license.
use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
use MOM_cpu_clock, only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE
use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
use MOM_diag_mediator, only : diag_ctrl, time_type
use MOM_debugging, only : hchksum, Bchksum
use MOM_error_handler, only : MOM_error, is_root_pe, FATAL, WARNING, NOTE
use MOM_file_parser, only : get_param, log_version, param_file_type
use MOM_grid, only : ocean_grid_type
use MOM_interface_heights, only : thickness_to_dz
use MOM_unit_scaling, only : unit_scale_type
use MOM_variables, only : thermo_var_ptrs
use MOM_verticalGrid, only : verticalGrid_type
use MOM_EOS, only : calculate_density_derivs
use MOM_EOS, only : calculate_density, calculate_specific_vol_derivs
implicit none ; private
#include <MOM_memory.h>
public Calculate_kappa_shear, Calc_kappa_shear_vertex, kappa_shear_init
public kappa_shear_is_used, kappa_shear_at_vertex
! 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.
!> This control structure holds the parameters that regulate shear mixing
type, public :: Kappa_shear_CS ; private
real :: RiNo_crit !< The critical shear Richardson number for
!! shear-entrainment [nondim]. The theoretical value is 0.25.
!! The values found by Jackson et al. are 0.25-0.35.
real :: Shearmix_rate !< A nondimensional rate scale for shear-driven
!! entrainment [nondim]. The value given by Jackson et al.
!! is 0.085-0.089.
real :: FRi_curvature !< A constant giving the curvature of the function
!! of the Richardson number that relates shear to
!! sources in the kappa equation [nondim].
!! The values found by Jackson et al. are -0.97 - -0.89.
real :: C_N !< The coefficient for the decay of TKE due to
!! stratification (i.e. proportional to N*tke) [nondim].
!! The values found by Jackson et al. are 0.24-0.28.
real :: C_S !< The coefficient for the decay of TKE due to
!! shear (i.e. proportional to |S|*tke) [nondim].
!! The values found by Jackson et al. are 0.14-0.12.
real :: lambda !< The coefficient for the buoyancy length scale
!! in the kappa equation [nondim].
!! The values found by Jackson et al. are 0.82-0.81.
real :: lambda2_N_S !< The square of the ratio of the coefficients of
!! the buoyancy and shear scales in the diffusivity
!! equation, 0 to eliminate the shear scale [nondim].
real :: lz_rescale !< A coefficient to rescale the distance to the nearest
!! solid boundary. This adjustment is to account for
!! regions where 3 dimensional turbulence prevents the
!! growth of shear instabilies [nondim].
real :: TKE_bg !< The background level of TKE [Z2 T-2 ~> m2 s-2].
real :: kappa_0 !< The background diapycnal diffusivity [H Z T-1 ~> m2 s-1 or Pa s]
real :: kappa_seed !< A moderately large seed value of diapycnal diffusivity that
!! is used as a starting turbulent diffusivity in the iterations
!! to findind an energetically constrained solution for the
!! shear-driven diffusivity [H Z T-1 ~> m2 s-1 or Pa s]
real :: kappa_trunc !< Diffusivities smaller than this are rounded to 0 [H Z T-1 ~> m2 s-1 or Pa s]
real :: kappa_tol_err !< The fractional error in kappa that is tolerated [nondim].
real :: Prandtl_turb !< Prandtl number used to convert Kd_shear into viscosity [nondim].
integer :: nkml !< The number of layers in the mixed layer, as
!! treated in this routine. If the pieces of the
!! mixed layer are not to be treated collectively,
!! nkml is set to 1.
integer :: max_RiNo_it !< The maximum number of iterations that may be used
!! to estimate the instantaneous shear-driven mixing.
integer :: max_KS_it !< The maximum number of iterations that may be used
!! to estimate the time-averaged diffusivity.
logical :: dKdQ_iteration_bug !< If true. use an older, dimensionally inconsistent estimate of
!! the derivative of diffusivity with energy in the Newton's method
!! iteration. The bug causes undercorrections when dz > 1m.
logical :: KS_at_vertex !< If true, do the calculations of the shear-driven mixing
!! at the cell vertices (i.e., the vorticity points).
logical :: eliminate_massless !< If true, massless layers are merged with neighboring
!! massive layers in this calculation.
! I can think of no good reason why this should be false. - RWH
real :: vel_underflow !< Velocity components smaller than vel_underflow
!! are set to 0 [L T-1 ~> m s-1].
real :: kappa_src_max_chg !< The maximum permitted increase in the kappa source within an
!! iteration relative to the local source [nondim]. This must be
!! greater than 1. The lower limit for the permitted fractional
!! decrease is (1 - 0.5/kappa_src_max_chg). These limits could
!! perhaps be made dynamic with an improved iterative solver.
logical :: psurf_bug !< If true, do a simple average of the cell surface pressures to get a
!! surface pressure at the corner if VERTEX_SHEAR=True. Otherwise mask
!! out any land points in the average.
logical :: all_layer_TKE_bug !< If true, report back the latest estimate of TKE instead of the
!! time average TKE when there is mass in all layers. Otherwise always
!! report the time-averaged TKE, as is currently done when there
!! are some massless layers.
logical :: restrictive_tolerance_check !< If false, uses the less restrictive tolerance check to
!! determine if a timestep is acceptable for the KS_it outer iteration
!! loop, as the code was originally written. True uses the more
!! restrictive check.
! logical :: layer_stagger = .false. ! If true, do the calculations centered at
! layers, rather than the interfaces.
logical :: debug = .false. !< If true, write verbose debugging messages.
type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to
!! regulate the timing of diagnostic output.
!>@{ Diagnostic IDs
integer :: id_Kd_shear = -1, id_TKE = -1
!>@}
end type Kappa_shear_CS
! integer :: id_clock_project, id_clock_KQ, id_clock_avg, id_clock_setup
contains
!> Subroutine for calculating shear-driven diffusivity and TKE in tracer columns
subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, &
kv_io, dt, G, GV, US, CS)
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
intent(in) :: u_in !< Initial zonal velocity [L T-1 ~> m s-1].
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
intent(in) :: v_in !< Initial meridional velocity [L T-1 ~> m s-1].
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any
!! available thermodynamic fields. Absent fields
!! have NULL ptrs.
real, dimension(:,:), pointer :: p_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa] (or NULL).
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
intent(inout) :: kappa_io !< The diapycnal diffusivity at each interface
!! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]. Initially this
!! is the value from the previous timestep, which may
!! accelerate the iteration toward convergence.
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
intent(out) :: tke_io !< The turbulent kinetic energy per unit mass at
!! each interface (not layer!) [Z2 T-2 ~> m2 s-2].
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
intent(inout) :: kv_io !< The vertical viscosity at each interface
!! (not layer!) [H Z T-1 ~> m2 s-1 or Pa s]. This discards any
!! previous value (i.e. it is intent out) and
!! simply sets Kv = Prandtl * Kd_shear
real, intent(in) :: dt !< Time increment [T ~> s].
type(Kappa_shear_CS), pointer :: CS !< The control structure returned by a previous
!! call to kappa_shear_init.
! Local variables
real, dimension(SZI_(G),SZK_(GV)) :: &
h_2d, & ! A 2-D version of h [H ~> m or kg m-2].
dz_2d, & ! Vertical distance between interface heights [Z ~> m].
u_2d, v_2d, & ! 2-D versions of u_in and v_in, converted to [L T-1 ~> m s-1].
T_2d, S_2d, rho_2d ! 2-D versions of T [C ~> degC], S [S ~> ppt], and rho [R ~> kg m-3].
real, dimension(SZI_(G),SZK_(GV)+1) :: &
kappa_2d, & ! 2-D version of kappa_io [H Z T-1 ~> m2 s-1 or Pa s]
tke_2d ! 2-D version tke_io [Z2 T-2 ~> m2 s-2].
real, dimension(SZK_(GV)) :: &
Idz, & ! The inverse of the thickness of the merged layers [H-1 ~> m2 kg-1].
h_lay, & ! The layer thickness [H ~> m or kg m-2]
dz_lay, & ! The geometric layer thickness in height units [Z ~> m]
u0xdz, & ! The initial zonal velocity times dz [H L T-1 ~> m2 s-1 or kg m-1 s-1]
v0xdz, & ! The initial meridional velocity times dz [H L T-1 ~> m2 s-1 or kg m-1 s-1]
T0xdz, & ! The initial temperature times thickness [C H ~> degC m or degC kg m-2] or if
! temperature is not a state variable, the density times thickness [R H ~> kg m-2 or kg2 m-3]
S0xdz ! The initial salinity times dz [S H ~> ppt m or ppt kg m-2].
real, dimension(SZK_(GV)+1) :: &
kappa, & ! The shear-driven diapycnal diffusivity at an interface [H Z T-1 ~> m2 s-1 or Pa s]
tke, & ! The Turbulent Kinetic Energy per unit mass at an interface [Z2 T-2 ~> m2 s-2].
kappa_avg, & ! The time-weighted average of kappa [H Z T-1 ~> m2 s-1 or Pa s]
tke_avg ! The time-weighted average of TKE [Z2 T-2 ~> m2 s-2].
real :: f2 ! The squared Coriolis parameter of each column [T-2 ~> s-2].
real :: surface_pres ! The top surface pressure [R L2 T-2 ~> Pa].
real :: dz_in_lay ! The running sum of the thickness in a layer [H ~> m or kg m-2]
real :: k0dt ! The background diffusivity times the timestep [H Z ~> m2 or kg m-1]
real :: dz_massless ! A layer thickness that is considered massless [H ~> m or kg m-2]
logical :: use_temperature ! If true, temperature and salinity have been
! allocated and are being used as state variables.
integer, dimension(SZK_(GV)+1) :: kc ! The index map between the original
! interfaces and the interfaces with massless layers
! merged into nearby massive layers.
real, dimension(SZK_(GV)+1) :: kf ! The fractional weight of interface kc+1 for
! interpolating back to the original index space [nondim].
integer :: is, ie, js, je, i, j, k, nz, nzc
is = G%isc ; ie = G%iec; js = G%jsc ; je = G%jec ; nz = GV%ke
use_temperature = associated(tv%T)
k0dt = dt*CS%kappa_0
dz_massless = 0.1*sqrt((US%Z_to_m*GV%m_to_H)*k0dt)
!$OMP parallel do default(private) shared(js,je,is,ie,nz,h,u_in,v_in,use_temperature,tv,G,GV,US, &
!$OMP CS,kappa_io,dz_massless,k0dt,p_surf,dt,tke_io,kv_io)
do j=js,je
! Convert layer thicknesses into geometric thickness in height units.
call thickness_to_dz(h, tv, dz_2d, j, G, GV)
do k=1,nz ; do i=is,ie
h_2d(i,k) = h(i,j,k)
u_2d(i,k) = u_in(i,j,k) ; v_2d(i,k) = v_in(i,j,k)
enddo ; enddo
if (use_temperature) then ; do k=1,nz ; do i=is,ie
T_2d(i,k) = tv%T(i,j,k) ; S_2d(i,k) = tv%S(i,j,k)
enddo ; enddo ; else ; do k=1,nz ; do i=is,ie
rho_2d(i,k) = GV%Rlay(k) ! Could be tv%Rho(i,j,k) ?
enddo ; enddo ; endif
!---------------------------------------
! Work on each column.
!---------------------------------------
do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then
! call cpu_clock_begin(id_clock_setup)
! Store a transposed version of the initial arrays.
! Any elimination of massless layers would occur here.
if (CS%eliminate_massless) then
nzc = 1
do k=1,nz
! Zero out the thicknesses of all layers, even if they are unused.
h_lay(k) = 0.0 ; dz_lay(k) = 0.0 ; u0xdz(k) = 0.0 ; v0xdz(k) = 0.0
T0xdz(k) = 0.0 ; S0xdz(k) = 0.0
! Add a new layer if this one has mass.
! if ((h_lay(nzc) > 0.0) .and. (h_2d(i,k) > dz_massless)) nzc = nzc+1
if ((k>CS%nkml) .and. (h_lay(nzc) > 0.0) .and. &
(h_2d(i,k) > dz_massless)) nzc = nzc+1
! Only merge clusters of massless layers.
! if ((h_lay(nzc) > dz_massless) .or. &
! ((h_lay(nzc) > 0.0) .and. (h_2d(i,k) > dz_massless))) nzc = nzc+1
kc(k) = nzc
h_lay(nzc) = h_lay(nzc) + h_2d(i,k)
dz_lay(nzc) = dz_lay(nzc) + dz_2d(i,k)
u0xdz(nzc) = u0xdz(nzc) + u_2d(i,k)*h_2d(i,k)
v0xdz(nzc) = v0xdz(nzc) + v_2d(i,k)*h_2d(i,k)
if (use_temperature) then
T0xdz(nzc) = T0xdz(nzc) + T_2d(i,k)*h_2d(i,k)
S0xdz(nzc) = S0xdz(nzc) + S_2d(i,k)*h_2d(i,k)
else
T0xdz(nzc) = T0xdz(nzc) + rho_2d(i,k)*h_2d(i,k)
S0xdz(nzc) = S0xdz(nzc) + rho_2d(i,k)*h_2d(i,k)
endif
enddo
kc(nz+1) = nzc+1
! Set up Idz as the inverse of layer thicknesses.
do k=1,nzc ; Idz(k) = 1.0 / h_lay(k) ; enddo
! Now determine kf, the fractional weight of interface kc when
! interpolating between interfaces kc and kc+1.
kf(1) = 0.0 ; dz_in_lay = h_2d(i,1)
do k=2,nz
if (kc(k) > kc(k-1)) then
kf(k) = 0.0 ; dz_in_lay = h_2d(i,k)
else
kf(k) = dz_in_lay*Idz(kc(k)) ; dz_in_lay = dz_in_lay + h_2d(i,k)
endif
enddo
kf(nz+1) = 0.0
else
do k=1,nz
h_lay(k) = h_2d(i,k)
dz_lay(k) = dz_2d(i,k)
u0xdz(k) = u_2d(i,k)*h_lay(k) ; v0xdz(k) = v_2d(i,k)*h_lay(k)
enddo
if (use_temperature) then
do k=1,nz
T0xdz(k) = T_2d(i,k)*h_lay(k) ; S0xdz(k) = S_2d(i,k)*h_lay(k)
enddo
else
do k=1,nz
T0xdz(k) = rho_2d(i,k)*h_lay(k) ; S0xdz(k) = rho_2d(i,k)*h_lay(k)
enddo
endif
nzc = nz
do k=1,nzc+1 ; kc(k) = k ; kf(k) = 0.0 ; enddo
endif
f2 = 0.25 * ((G%CoriolisBu(I,j)**2 + G%CoriolisBu(I-1,J-1)**2) + &
(G%CoriolisBu(I,J-1)**2 + G%CoriolisBu(I-1,J)**2))
surface_pres = 0.0 ; if (associated(p_surf)) surface_pres = p_surf(i,j)
! ---------------------------------------------------- I_Ld2_1d, dz_Int_1d
! Set the initial guess for kappa, here defined at interfaces.
! ----------------------------------------------------
do K=1,nzc+1 ; kappa(K) = CS%kappa_seed ; enddo
call kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, &
h_lay, dz_lay, u0xdz, v0xdz, T0xdz, S0xdz, kappa_avg, &
tke_avg, tv, CS, GV, US)
! call cpu_clock_begin(id_clock_setup)
! Extrapolate from the vertically reduced grid back to the original layers.
if (nz == nzc) then
do K=1,nz+1
kappa_2d(i,K) = kappa_avg(K)
if (CS%all_layer_TKE_bug) then
tke_2d(i,K) = tke(K)
else
tke_2d(i,K) = tke_avg(K)
endif
enddo
else
do K=1,nz+1
if (kf(K) == 0.0) then
kappa_2d(i,K) = kappa_avg(kc(K))
tke_2d(i,K) = tke_avg(kc(K))
else
kappa_2d(i,K) = (1.0-kf(K)) * kappa_avg(kc(K)) + &
kf(K) * kappa_avg(kc(K)+1)
tke_2d(i,K) = (1.0-kf(K)) * tke_avg(kc(K)) + &
kf(K) * tke_avg(kc(K)+1)
endif
enddo
endif
! call cpu_clock_end(id_clock_setup)
else ! Land points, still inside the i-loop.
do K=1,nz+1
kappa_2d(i,K) = 0.0 ; tke_2d(i,K) = 0.0
enddo
endif ; enddo ! i-loop
do K=1,nz+1 ; do i=is,ie
kappa_io(i,j,K) = G%mask2dT(i,j) * kappa_2d(i,K)
tke_io(i,j,K) = G%mask2dT(i,j) * tke_2d(i,K)
kv_io(i,j,K) = ( G%mask2dT(i,j) * kappa_2d(i,K) ) * CS%Prandtl_turb
enddo ; enddo
enddo ! end of j-loop
if (CS%debug) then
call hchksum(kappa_io, "kappa", G%HI, unscale=GV%HZ_T_to_m2_s)
call hchksum(tke_io, "tke", G%HI, unscale=US%Z_to_m**2*US%s_to_T**2)
endif
if (CS%id_Kd_shear > 0) call post_data(CS%id_Kd_shear, kappa_io, CS%diag)
if (CS%id_TKE > 0) call post_data(CS%id_TKE, tke_io, CS%diag)
end subroutine Calculate_kappa_shear
!> Subroutine for calculating shear-driven diffusivity and TKE in corner columns
subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_io, tke_io, &
kv_io, dt, G, GV, US, CS)
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
type(verticalGrid_type), intent(in) :: GV !< The ocean's 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_in !< Initial zonal velocity [L T-1 ~> m s-1].
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
intent(in) :: v_in !< Initial meridional velocity [L T-1 ~> m s-1].
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)), &
intent(in) :: T_in !< Layer potential temperatures [C ~> degC]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
intent(in) :: S_in !< Layer salinities [S ~> ppt]
type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any
!! available thermodynamic fields. Absent fields
!! have NULL ptrs.
real, dimension(:,:), pointer :: p_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa]
!! (or NULL).
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
intent(out) :: kappa_io !< The diapycnal diffusivity at each interface
!! (not layer!) [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)+1), &
intent(out) :: tke_io !< The turbulent kinetic energy per unit mass at
!! each interface (not layer!) [Z2 T-2 ~> m2 s-2].
real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)+1), &
intent(inout) :: kv_io !< The vertical viscosity at each interface
!! [H Z T-1 ~> m2 s-1 or Pa s].
!! The previous value is used to initialize kappa
!! in the vertex columns as Kappa = Kv/Prandtl
!! to accelerate the iteration toward convergence.
real, intent(in) :: dt !< Time increment [T ~> s].
type(Kappa_shear_CS), pointer :: CS !< The control structure returned by a previous
!! call to kappa_shear_init.
! Local variables
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
dz_3d ! Vertical distance between interface heights [Z ~> m].
real, dimension(SZIB_(G),SZK_(GV)) :: &
h_2d, & ! A 2-D version of h [H ~> m or kg m-2].
dz_2d, & ! Vertical distance between interface heights [Z ~> m].
u_2d, v_2d, & ! 2-D versions of u_in and v_in, converted to [L T-1 ~> m s-1].
T_2d, S_2d, rho_2d ! 2-D versions of T [C ~> degC], S [S ~> ppt], and rho [R ~> kg m-3].
real, dimension(SZIB_(G),SZK_(GV)+1,2) :: &
kappa_2d ! Quasi 2-D versions of kappa_io [H Z T-1 ~> m2 s-1 or Pa s]
real, dimension(SZIB_(G),SZK_(GV)+1) :: &
tke_2d ! 2-D version tke_io [Z2 T-2 ~> m2 s-2].
real, dimension(SZK_(GV)) :: &
Idz, & ! The inverse of the thickness of the merged layers [H-1 ~> m2 kg-1].
h_lay, & ! The layer thickness [H ~> m or kg m-2]
dz_lay, & ! The geometric layer thickness in height units [Z ~> m]
u0xdz, & ! The initial zonal velocity times dz [L H T-1 ~> m2 s-1 or kg m-1 s-1].
v0xdz, & ! The initial meridional velocity times dz [H L T-1 ~> m2 s-1 or kg m-1 s-1]
T0xdz, & ! The initial temperature times dz [C H ~> degC m or degC kg m-2]
S0xdz ! The initial salinity times dz [S H ~> ppt m or ppt kg m-2]
real, dimension(SZK_(GV)+1) :: &
kappa, & ! The shear-driven diapycnal diffusivity at an interface [H Z T-1 ~> m2 s-1 or Pa s]
tke, & ! The Turbulent Kinetic Energy per unit mass at an interface [Z2 T-2 ~> m2 s-2].
kappa_avg, & ! The time-weighted average of kappa [H Z T-1 ~> m2 s-1 or Pa s]
tke_avg ! The time-weighted average of TKE [Z2 T-2 ~> m2 s-2].
real :: f2 ! The squared Coriolis parameter of each column [T-2 ~> s-2].
real :: surface_pres ! The top surface pressure [R L2 T-2 ~> Pa].
real :: dz_in_lay ! The running sum of the thickness in a layer [H ~> m or kg m-2]
real :: k0dt ! The background diffusivity times the timestep [H Z ~> m2 or kg m-1]
real :: dz_massless ! A layer thickness that is considered massless [H ~> m or kg m-2]
real :: I_hwt ! The inverse of the masked thickness weights [H-1 ~> m-1 or m2 kg-1].
real :: I_Prandtl ! The inverse of the turbulent Prandtl number [nondim].
logical :: use_temperature ! If true, temperature and salinity have been
! allocated and are being used as state variables.
integer, dimension(SZK_(GV)+1) :: kc ! The index map between the original
! interfaces and the interfaces with massless layers
! merged into nearby massive layers.
real, dimension(SZK_(GV)+1) :: kf ! The fractional weight of interface kc+1 for
! interpolating back to the original index space [nondim].
integer :: IsB, IeB, JsB, JeB, i, j, k, nz, nzc, J2, J2m1
! Diagnostics that should be deleted?
isB = G%isc-1 ; ieB = G%iecB ; jsB = G%jsc-1 ; jeB = G%jecB ; nz = GV%ke
use_temperature = associated(tv%T)
k0dt = dt*CS%kappa_0
dz_massless = 0.1*sqrt((US%Z_to_m*GV%m_to_H)*k0dt)
I_Prandtl = 0.0 ; if (CS%Prandtl_turb > 0.0) I_Prandtl = 1.0 / CS%Prandtl_turb
! Convert layer thicknesses into geometric thickness in height units.
call thickness_to_dz(h, tv, dz_3d, G, GV, US, halo_size=1)
!$OMP parallel do default(private) shared(jsB,jeB,isB,ieB,nz,h,u_in,v_in,use_temperature,tv,G,GV, &
!$OMP US,CS,kappa_io,dz_massless,k0dt,p_surf,dt,tke_io,kv_io,I_Prandtl)
do J=JsB,JeB
J2 = mod(J,2)+1 ; J2m1 = 3-J2 ! = mod(J-1,2)+1
! Interpolate the various quantities to the corners, using masks.
do k=1,nz ; do I=IsB,IeB
u_2d(I,k) = (u_in(I,j,k) * (G%mask2dCu(I,j) * (h(i,j,k) + h(i+1,j,k))) + &
u_in(I,j+1,k) * (G%mask2dCu(I,j+1) * (h(i,j+1,k) + h(i+1,j+1,k))) ) / &
((G%mask2dCu(I,j) * (h(i,j,k) + h(i+1,j,k)) + &
G%mask2dCu(I,j+1) * (h(i,j+1,k) + h(i+1,j+1,k))) + GV%H_subroundoff)
v_2d(I,k) = (v_in(i,J,k) * (G%mask2dCv(i,J) * (h(i,j,k) + h(i,j+1,k))) + &
v_in(i+1,J,k) * (G%mask2dCv(i+1,J) * (h(i+1,j,k) + h(i+1,j+1,k))) ) / &
((G%mask2dCv(i,J) * (h(i,j,k) + h(i,j+1,k)) + &
G%mask2dCv(i+1,J) * (h(i+1,j,k) + h(i+1,j+1,k))) + GV%H_subroundoff)
I_hwt = 1.0 / (((G%mask2dT(i,j) * h(i,j,k) + G%mask2dT(i+1,j+1) * h(i+1,j+1,k)) + &
(G%mask2dT(i+1,j) * h(i+1,j,k) + G%mask2dT(i,j+1) * h(i,j+1,k))) + &
GV%H_subroundoff)
if (use_temperature) then
T_2d(I,k) = ( ((G%mask2dT(i,j) * h(i,j,k)) * T_in(i,j,k) + &
(G%mask2dT(i+1,j+1) * h(i+1,j+1,k)) * T_in(i+1,j+1,k)) + &
((G%mask2dT(i+1,j) * h(i+1,j,k)) * T_in(i+1,j,k) + &
(G%mask2dT(i,j+1) * h(i,j+1,k)) * T_in(i,j+1,k)) ) * I_hwt
S_2d(I,k) = ( ((G%mask2dT(i,j) * h(i,j,k)) * S_in(i,j,k) + &
(G%mask2dT(i+1,j+1) * h(i+1,j+1,k)) * S_in(i+1,j+1,k)) + &
((G%mask2dT(i+1,j) * h(i+1,j,k)) * S_in(i+1,j,k) + &
(G%mask2dT(i,j+1) * h(i,j+1,k)) * S_in(i,j+1,k)) ) * I_hwt
endif
h_2d(I,k) = ((G%mask2dT(i,j) * h(i,j,k) + G%mask2dT(i+1,j+1) * h(i+1,j+1,k)) + &
(G%mask2dT(i+1,j) * h(i+1,j,k) + G%mask2dT(i,j+1) * h(i,j+1,k)) ) / &
((G%mask2dT(i,j) + G%mask2dT(i+1,j+1)) + &
(G%mask2dT(i+1,j) + G%mask2dT(i,j+1)) + 1.0e-36 )
dz_2d(I,k) = ((G%mask2dT(i,j) * dz_3d(i,j,k) + G%mask2dT(i+1,j+1) * dz_3d(i+1,j+1,k)) + &
(G%mask2dT(i+1,j) * dz_3d(i+1,j,k) + G%mask2dT(i,j+1) * dz_3d(i,j+1,k)) ) / &
((G%mask2dT(i,j) + G%mask2dT(i+1,j+1)) + &
(G%mask2dT(i+1,j) + G%mask2dT(i,j+1)) + 1.0e-36 )
! h_2d(I,k) = 0.25*((h(i,j,k) + h(i+1,j+1,k)) + (h(i+1,j,k) + h(i,j+1,k)))
! h_2d(I,k) = ((h(i,j,k)**2 + h(i+1,j+1,k)**2) + &
! (h(i+1,j,k)**2 + h(i,j+1,k)**2)) * I_hwt
enddo ; enddo
if (.not.use_temperature) then ; do k=1,nz ; do I=IsB,IeB
rho_2d(I,k) = GV%Rlay(k)
enddo ; enddo ; endif
!---------------------------------------
! Work on each column.
!---------------------------------------
do I=IsB,IeB ; if ((G%mask2dCu(I,j) + G%mask2dCu(I,j+1)) + &
(G%mask2dCv(i,J) + G%mask2dCv(i+1,J)) > 0.0) then
! call cpu_clock_begin(Id_clock_setup)
! Store a transposed version of the initial arrays.
! Any elimination of massless layers would occur here.
if (CS%eliminate_massless) then
nzc = 1
do k=1,nz
! Zero out the thicknesses of all layers, even if they are unused.
h_lay(k) = 0.0 ; dz_lay(k) = 0.0 ; u0xdz(k) = 0.0 ; v0xdz(k) = 0.0
T0xdz(k) = 0.0 ; S0xdz(k) = 0.0
! Add a new layer if this one has mass.
! if ((h_lay(nzc) > 0.0) .and. (h_2d(I,k) > dz_massless)) nzc = nzc+1
if ((k>CS%nkml) .and. (h_lay(nzc) > 0.0) .and. &
(h_2d(I,k) > dz_massless)) nzc = nzc+1
! Only merge clusters of massless layers.
! if ((h_lay(nzc) > dz_massless) .or. &
! ((h_lay(nzc) > 0.0) .and. (h_2d(I,k) > dz_massless))) nzc = nzc+1
kc(k) = nzc
h_lay(nzc) = h_lay(nzc) + h_2d(I,k)
dz_lay(nzc) = dz_lay(nzc) + dz_2d(I,k)
u0xdz(nzc) = u0xdz(nzc) + u_2d(I,k)*h_2d(I,k)
v0xdz(nzc) = v0xdz(nzc) + v_2d(I,k)*h_2d(I,k)
if (use_temperature) then
T0xdz(nzc) = T0xdz(nzc) + T_2d(I,k)*h_2d(I,k)
S0xdz(nzc) = S0xdz(nzc) + S_2d(I,k)*h_2d(I,k)
else
T0xdz(nzc) = T0xdz(nzc) + rho_2d(I,k)*h_2d(I,k)
S0xdz(nzc) = S0xdz(nzc) + rho_2d(I,k)*h_2d(I,k)
endif
enddo
kc(nz+1) = nzc+1
! Set up Idz as the inverse of layer thicknesses.
do k=1,nzc ; Idz(k) = 1.0 / h_lay(k) ; enddo
! Now determine kf, the fractional weight of interface kc when
! interpolating between interfaces kc and kc+1.
kf(1) = 0.0 ; dz_in_lay = h_2d(I,1)
do k=2,nz
if (kc(k) > kc(k-1)) then
kf(k) = 0.0 ; dz_in_lay = h_2d(I,k)
else
kf(k) = dz_in_lay*Idz(kc(k)) ; dz_in_lay = dz_in_lay + h_2d(I,k)
endif
enddo
kf(nz+1) = 0.0
else
do k=1,nz
h_lay(k) = h_2d(I,k)
dz_lay(k) = dz_2d(I,k)
u0xdz(k) = u_2d(I,k)*h_lay(k) ; v0xdz(k) = v_2d(I,k)*h_lay(k)
enddo
if (use_temperature) then
do k=1,nz
T0xdz(k) = T_2d(I,k)*h_lay(k) ; S0xdz(k) = S_2d(I,k)*h_lay(k)
enddo
else
do k=1,nz
T0xdz(k) = rho_2d(I,k)*h_lay(k) ; S0xdz(k) = rho_2d(I,k)*h_lay(k)
enddo
endif
nzc = nz
do k=1,nzc+1 ; kc(k) = k ; kf(k) = 0.0 ; enddo
endif
f2 = G%CoriolisBu(I,J)**2
surface_pres = 0.0
if (associated(p_surf)) then
if (CS%psurf_bug) then
! This is wrong because it is averaging values from land in some places.
surface_pres = 0.25 * ((p_surf(i,j) + p_surf(i+1,j+1)) + &
(p_surf(i+1,j) + p_surf(i,j+1)))
else
surface_pres = ((G%mask2dT(i,j) * p_surf(i,j) + G%mask2dT(i+1,j+1) * p_surf(i+1,j+1)) + &
(G%mask2dT(i+1,j) * p_surf(i+1,j) + G%mask2dT(i,j+1) * p_surf(i,j+1)) ) / &
((G%mask2dT(i,j) + G%mask2dT(i+1,j+1)) + &
(G%mask2dT(i+1,j) + G%mask2dT(i,j+1)) + 1.0e-36 )
endif
endif
! ----------------------------------------------------
! Set the initial guess for kappa, here defined at interfaces.
! ----------------------------------------------------
do K=1,nzc+1 ; kappa(K) = CS%kappa_seed ; enddo
call kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, &
h_lay, dz_lay, u0xdz, v0xdz, T0xdz, S0xdz, kappa_avg, &
tke_avg, tv, CS, GV, US)
! call cpu_clock_begin(Id_clock_setup)
! Extrapolate from the vertically reduced grid back to the original layers.
if (nz == nzc) then
do K=1,nz+1
kappa_2d(I,K,J2) = kappa_avg(K)
if (CS%all_layer_TKE_bug) then
tke_2d(i,K) = tke(K)
else
tke_2d(i,K) = tke_avg(K)
endif
enddo
else
do K=1,nz+1
if (kf(K) == 0.0) then
kappa_2d(I,K,J2) = kappa_avg(kc(K))
tke_2d(I,K) = tke_avg(kc(K))
else
kappa_2d(I,K,J2) = (1.0-kf(K)) * kappa_avg(kc(K)) + kf(K) * kappa_avg(kc(K)+1)
tke_2d(I,K) = (1.0-kf(K)) * tke_avg(kc(K)) + kf(K) * tke_avg(kc(K)+1)
endif
enddo
endif
! call cpu_clock_end(Id_clock_setup)
else ! Land points, still inside the i-loop.
do K=1,nz+1
kappa_2d(I,K,J2) = 0.0 ; tke_2d(I,K) = 0.0
enddo
endif ; enddo ! i-loop
do K=1,nz+1 ; do I=IsB,IeB
tke_io(I,J,K) = G%mask2dBu(I,J) * tke_2d(I,K)
kv_io(I,J,K) = ( G%mask2dBu(I,J) * kappa_2d(I,K,J2) ) * CS%Prandtl_turb
enddo ; enddo
if (J>=G%jsc) then ; do K=1,nz+1 ; do i=G%isc,G%iec
! Set the diffusivities in tracer columns from the values at vertices.
kappa_io(i,j,K) = G%mask2dT(i,j) * 0.25 * &
((kappa_2d(I-1,K,J2m1) + kappa_2d(I,K,J2)) + &
(kappa_2d(I-1,K,J2) + kappa_2d(I,K,J2m1)))
enddo ; enddo ; endif
enddo ! end of J-loop
if (CS%debug) then
call hchksum(kappa_io, "kappa", G%HI, unscale=GV%HZ_T_to_m2_s)
call Bchksum(tke_io, "tke", G%HI, unscale=US%Z_to_m**2*US%s_to_T**2)
endif
if (CS%id_Kd_shear > 0) call post_data(CS%id_Kd_shear, kappa_io, CS%diag)
if (CS%id_TKE > 0) call post_data(CS%id_TKE, tke_io, CS%diag)
end subroutine Calc_kappa_shear_vertex
!> This subroutine calculates shear-driven diffusivity and TKE in a single column
subroutine kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, hlay, dz_lay, &
u0xdz, v0xdz, T0xdz, S0xdz, kappa_avg, tke_avg, tv, CS, GV, US)
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
real, dimension(SZK_(GV)+1), &
intent(inout) :: kappa !< The time-weighted average of kappa [H Z T-1 ~> m2 s-1 or Pa s]
real, dimension(SZK_(GV)+1), &
intent(out) :: tke !< The Turbulent Kinetic Energy per unit mass at
!! an interface [Z2 T-2 ~> m2 s-2].
integer, intent(in) :: nzc !< The number of active layers in the column.
real, intent(in) :: f2 !< The square of the Coriolis parameter [T-2 ~> s-2].
real, intent(in) :: surface_pres !< The surface pressure [R L2 T-2 ~> Pa].
real, dimension(SZK_(GV)), &
intent(in) :: hlay !< The layer thickness [H ~> m or kg m-2]
real, dimension(SZK_(GV)), &
intent(in) :: dz_lay !< The geometric layer thickness in height units [Z ~> m]
real, dimension(SZK_(GV)), &
intent(in) :: u0xdz !< The initial zonal velocity times hlay [H L T-1 ~> m2 s-1 or kg m-1 s-1]
real, dimension(SZK_(GV)), &
intent(in) :: v0xdz !< The initial meridional velocity times the
!! layer thickness [H L T-1 ~> m2 s-1 or kg m-1 s-1]
real, dimension(SZK_(GV)), &
intent(in) :: T0xdz !< The initial temperature times hlay [C H ~> degC m or degC kg m-2]
real, dimension(SZK_(GV)), &
intent(in) :: S0xdz !< The initial salinity times hlay [S H ~> ppt m or ppt kg m-2]
real, dimension(SZK_(GV)+1), &
intent(out) :: kappa_avg !< The time-weighted average of kappa [H Z T-1 ~> m2 s-1 or Pa s]
real, dimension(SZK_(GV)+1), &
intent(out) :: tke_avg !< The time-weighted average of TKE [Z2 T-2 ~> m2 s-2].
real, intent(in) :: dt !< Time increment [T ~> s].
type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any
!! available thermodynamic fields. Absent fields
!! have NULL ptrs.
type(Kappa_shear_CS), pointer :: CS !< The control structure returned by a previous
!! call to kappa_shear_init.
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
! Local variables
real, dimension(nzc) :: &
u, & ! The zonal velocity after a timestep of mixing [L T-1 ~> m s-1].
v, & ! The meridional velocity after a timestep of mixing [L T-1 ~> m s-1].
Idz, & ! The inverse of the distance between TKE points [Z-1 ~> m-1].
T, & ! The potential temperature after a timestep of mixing [C ~> degC].
Sal, & ! The salinity after a timestep of mixing [S ~> ppt].
u_test, v_test, & ! Temporary velocities [L T-1 ~> m s-1].
T_test, S_test ! Temporary temperatures [C ~> degC] and salinities [S ~> ppt].
real, dimension(nzc+1) :: &
N2, & ! The squared buoyancy frequency at an interface [T-2 ~> s-2].
h_Int, & ! The extent of a finite-volume space surrounding an interface,
! as used in calculating kappa and TKE [H ~> m or kg m-2]
dz_Int, & ! The vertical distance with the space surrounding an interface,
! as used in calculating kappa and TKE [Z ~> m]
dz_h_Int, & ! The ratio of the vertical distances to the thickness around an
! interface [Z H-1 ~> nondim or m3 kg-1]. In non-Boussinesq mode
! this is the specific volume, otherwise it is a scaling factor.
I_dz_int, & ! The inverse of the distance between velocity & density points
! above and below an interface [Z-1 ~> m-1]. This is used to
! calculate N2, shear and fluxes.
S2, & ! The squared shear at an interface [T-2 ~> s-2].
a1, & ! a1 is the coupling between adjacent interfaces in the TKE,
! velocity, and density equations [H ~> m or kg m-2]
c1, & ! c1 is used in the tridiagonal (and similar) solvers [nondim].
k_src, & ! The shear-dependent source term in the kappa equation [T-1 ~> s-1]
kappa_src, & ! The shear-dependent source term in the kappa equation [T-1 ~> s-1]
kappa_out, & ! The kappa that results from the kappa equation [H Z T-1 ~> m2 s-1 or Pa s]
kappa_mid, & ! The average of the initial and predictor estimates of kappa [H Z T-1 ~> m2 s-1 or Pa s]
tke_pred, & ! The value of TKE from a predictor step [Z2 T-2 ~> m2 s-2].
kappa_pred, & ! The value of kappa from a predictor step [H Z T-1 ~> m2 s-1 or Pa s]
pressure, & ! The pressure at an interface [R L2 T-2 ~> Pa].
T_int, & ! The temperature interpolated to an interface [C ~> degC].
Sal_int, & ! The salinity interpolated to an interface [S ~> ppt].
dbuoy_dT, & ! The partial derivative of buoyancy with changes in temperature [Z T-2 C-1 ~> m s-2 degC-1]
dbuoy_dS, & ! The partial derivative of buoyancy with changes in salinity [Z T-2 S-1 ~> m s-2 ppt-1]
dSpV_dT, & ! The partial derivative of specific volume with changes in temperature [R-1 C-1 ~> m3 kg-1 degC-1]
dSpV_dS, & ! The partial derivative of specific volume with changes in salinity [R-1 S-1 ~> m3 kg-1 ppt-1]
rho_int, & ! The in situ density interpolated to an interface [R ~> kg m-3]
I_L2_bdry, & ! The inverse of the square of twice the harmonic mean
! distance to the top and bottom boundaries [H-1 Z-1 ~> m-2 or m kg-1].
K_Q, & ! Diffusivity divided by TKE [H T Z-1 ~> s or kg s m-3]
K_Q_tmp, & ! A temporary copy of diffusivity divided by TKE [H T Z-1 ~> s or kg s m-3]
local_src_avg, & ! The time-integral of the local source [nondim]
tol_min, & ! Minimum tolerated ksrc for the corrector step [T-1 ~> s-1].
tol_max, & ! Maximum tolerated ksrc for the corrector step [T-1 ~> s-1].
tol_chg, & ! The tolerated kappa change integrated over a timestep [nondim].
dist_from_top, & ! The distance from the top surface [Z ~> m].
h_from_top, & ! The total thickness above an interface [H ~> m or kg m-2]
local_src ! The sum of all sources of kappa, including kappa_src and
! sources from the elliptic term [T-1 ~> s-1]
real :: dist_from_bot ! The distance from the bottom surface [Z ~> m].
real :: h_from_bot ! The total thickness below and interface [H ~> m or kg m-2]
real :: b1 ! The inverse of the pivot in the tridiagonal equations [H-1 ~> m-1 or m2 kg-1].
real :: bd1 ! A term in the denominator of b1 [H ~> m or kg m-2].
real :: d1 ! 1 - c1 in the tridiagonal equations [nondim]
real :: gR0 ! A conversion factor from H to pressure, Rho_0 times g in Boussinesq
! mode, or just g when non-Boussinesq [R L2 T-2 H-1 ~> kg m-2 s-2 or m s-2].
real :: g_R0 ! g_R0 is a rescaled version of g/Rho [Z R-1 T-2 ~> m4 kg-1 s-2].
real :: Norm ! A factor that normalizes two weights to 1 [H-2 ~> m-2 or m4 kg-2].
real :: tol_dksrc ! Tolerance for the change in the kappa source within an iteration
! relative to the local source [nondim]. This must be greater than 1.
real :: tol2 ! The tolerance for the change in the kappa source within an iteration
! relative to the average local source over previous iterations [nondim].
real :: tol_dksrc_low ! The tolerance for the fractional decrease in ksrc
! within an iteration [nondim]. 0 < tol_dksrc_low < 1.
real :: Ri_crit ! The critical shear Richardson number for shear-
! driven mixing [nondim]. The theoretical value is 0.25.
real :: dt_rem ! The remaining time to advance the solution [T ~> s].
real :: dt_now ! The time step used in the current iteration [T ~> s].
real :: dt_wt ! The fractional weight of the current iteration [nondim].
real :: dt_test ! A time-step that is being tested for whether it
! gives acceptably small changes in k_src [T ~> s].
real :: Idtt ! Idtt = 1 / dt_test [T-1 ~> s-1].
real :: dt_inc ! An increment to dt_test that is being tested [T ~> s].
real :: wt_a ! The fraction of a layer thickness identified with the interface
! above a layer [nondim]
real :: wt_b ! The fraction of a layer thickness identified with the interface
! below a layer [nondim]
real :: k0dt ! The background diffusivity times the timestep [H Z ~> m2 or kg m-1].
real :: I_lz_rescale_sqr ! The inverse of a rescaling factor for L2_bdry (Lz) squared [nondim].
logical :: valid_dt ! If true, all levels so far exhibit acceptably small changes in k_src.
logical :: use_temperature ! If true, temperature and salinity have been
! allocated and are being used as state variables.
integer :: ks_kappa, ke_kappa ! The k-range with nonzero kappas.
integer :: dt_refinements ! The number of 2-fold refinements that will be used
! to estimate the maximum permitted time step. I.e.,
! the resolution is 1/2^dt_refinements.
integer :: k, itt, itt_dt
! This calculation of N2 is for debugging only.
! real, dimension(SZK_(GV)+1) :: &
! N2_debug, & ! A version of N2 for debugging [T-2 ~> s-2]
Ri_crit = CS%Rino_crit
gR0 = GV%H_to_RZ * GV%g_Earth
g_R0 = (US%L_to_Z**2 * GV%g_Earth) / GV%Rho0
k0dt = dt*CS%kappa_0
I_lz_rescale_sqr = 1.0; if (CS%lz_rescale > 0) I_lz_rescale_sqr = 1/(CS%lz_rescale*CS%lz_rescale)
tol_dksrc = CS%kappa_src_max_chg
if (tol_dksrc == 10.0) then
! This is equivalent to the expression below, but avoids changes at roundoff for the default value.
tol_dksrc_low = 0.95
else
tol_dksrc_low = (tol_dksrc - 0.5)/tol_dksrc
endif
tol2 = 2.0*CS%kappa_tol_err
dt_refinements = 5 ! Selected so that 1/2^dt_refinements < 1-tol_dksrc_low
use_temperature = .false. ; if (associated(tv%T)) use_temperature = .true.
! Set up Idz as the inverse of layer thicknesses.
do k=1,nzc ; Idz(k) = 1.0 / dz_lay(k) ; enddo
! Set up I_dz_int as the inverse of the distance between
! adjacent layer centers.
I_dz_int(1) = 2.0 / dz_lay(1)
dist_from_top(1) = 0.0 ; h_from_top(1) = 0.0
do K=2,nzc
I_dz_int(K) = 2.0 / (dz_lay(k-1) + dz_lay(k))
dist_from_top(K) = dist_from_top(K-1) + dz_lay(k-1)
h_from_top(K) = h_from_top(K-1) + hlay(k-1)
enddo
I_dz_int(nzc+1) = 2.0 / dz_lay(nzc)
! Find the inverse of the squared distances from the boundaries.
dist_from_bot = 0.0 ; h_from_bot = 0.0
do K=nzc,2,-1
dist_from_bot = dist_from_bot + dz_lay(k)
h_from_bot = h_from_bot + hlay(k)
! Find the inverse of the squared distances from the boundaries,
I_L2_bdry(K) = ((dist_from_top(K) + dist_from_bot) * (h_from_top(K) + h_from_bot)) / &
((dist_from_top(K) * dist_from_bot) * (h_from_top(K) * h_from_bot))
! reduce the distance by a factor of "lz_rescale"
I_L2_bdry(K) = I_lz_rescale_sqr*I_L2_bdry(K)
enddo
! Determine the velocities and thicknesses after eliminating massless
! layers and applying a time-step of background diffusion.
if (nzc > 1) then
a1(2) = k0dt*I_dz_int(2)
b1 = 1.0 / (hlay(1) + a1(2))
u(1) = b1 * u0xdz(1) ; v(1) = b1 * v0xdz(1)
T(1) = b1 * T0xdz(1) ; Sal(1) = b1 * S0xdz(1)
c1(2) = a1(2) * b1 ; d1 = hlay(1) * b1 ! = 1 - c1
do k=2,nzc-1
bd1 = hlay(k) + d1*a1(k)
a1(k+1) = k0dt*I_dz_int(k+1)
b1 = 1.0 / (bd1 + a1(k+1))
u(k) = b1 * (u0xdz(k) + a1(k)*u(k-1))
v(k) = b1 * (v0xdz(k) + a1(k)*v(k-1))
T(k) = b1 * (T0xdz(k) + a1(k)*T(k-1))
Sal(k) = b1 * (S0xdz(k) + a1(k)*Sal(k-1))
c1(k+1) = a1(k+1) * b1 ; d1 = bd1 * b1 ! d1 = 1 - c1
enddo
! rho or T and S have insulating boundary conditions, u & v use no-slip
! bottom boundary conditions (if kappa0 > 0).
! For no-slip bottom boundary conditions
b1 = 1.0 / ((hlay(nzc) + d1*a1(nzc)) + k0dt*I_dz_int(nzc+1))
u(nzc) = b1 * (u0xdz(nzc) + a1(nzc)*u(nzc-1))
v(nzc) = b1 * (v0xdz(nzc) + a1(nzc)*v(nzc-1))
! For insulating boundary conditions
b1 = 1.0 / (hlay(nzc) + d1*a1(nzc))
T(nzc) = b1 * (T0xdz(nzc) + a1(nzc)*T(nzc-1))
Sal(nzc) = b1 * (S0xdz(nzc) + a1(nzc)*Sal(nzc-1))
do k=nzc-1,1,-1
u(k) = u(k) + c1(k+1)*u(k+1) ; v(k) = v(k) + c1(k+1)*v(k+1)
T(k) = T(k) + c1(k+1)*T(k+1) ; Sal(k) = Sal(k) + c1(k+1)*Sal(k+1)
enddo
else
! This is correct, but probably unnecessary.
b1 = 1.0 / (hlay(1) + k0dt*I_dz_int(2))
u(1) = b1 * u0xdz(1) ; v(1) = b1 * v0xdz(1)
b1 = 1.0 / hlay(1)
T(1) = b1 * T0xdz(1) ; Sal(1) = b1 * S0xdz(1)
endif
! This uses half the harmonic mean of thicknesses to provide two estimates
! of the boundary between cells, and the inverse of the harmonic mean to
! weight the two estimates. The net effect is that interfaces around thin
! layers have thin cells, and the total thickness adds up properly.
! The top- and bottom- interfaces have zero thickness, consistent with
! adding additional zero thickness layers.
h_Int(1) = 0.0 ; h_Int(2) = hlay(1)
dz_Int(1) = 0.0 ; dz_Int(2) = dz_lay(1)
do K=2,nzc-1
Norm = 1.0 / (hlay(k)*(hlay(k-1)+hlay(k+1)) + 2.0*hlay(k-1)*hlay(k+1))
wt_a = ((hlay(k)+hlay(k+1)) * hlay(k-1)) * Norm
wt_b = ((hlay(k-1)+hlay(k)) * hlay(k+1)) * Norm
h_Int(K) = h_Int(K) + hlay(k) * wt_a
h_Int(K+1) = hlay(k) * wt_b
dz_Int(K) = dz_Int(K) + dz_lay(k) * wt_a
dz_Int(K+1) = dz_lay(k) * wt_b
enddo
h_Int(nzc) = h_Int(nzc) + hlay(nzc) ; h_Int(nzc+1) = 0.0
dz_Int(nzc) = dz_Int(nzc) + dz_lay(nzc) ; dz_Int(nzc+1) = 0.0
if (GV%Boussinesq) then
do K=1,nzc+1 ; dz_h_Int(K) = GV%H_to_Z ; enddo
else
! Find an effective average specific volume around an interface.
dz_h_Int(1:nzc+1) = 0.0
if (hlay(1) > 0.0) dz_h_Int(1) = dz_lay(1) / hlay(1)
do K=2,nzc+1
if (h_Int(K) > 0.0) then
dz_h_Int(K) = dz_Int(K) / h_Int(K)
else
dz_h_Int(K) = dz_h_Int(K-1)
endif
enddo
endif
! Calculate thermodynamic coefficients and an initial estimate of N2.
if (use_temperature) then
pressure(1) = surface_pres
do K=2,nzc
pressure(K) = pressure(K-1) + gR0*hlay(k-1)
T_int(K) = 0.5*(T(k-1) + T(k))
Sal_int(K) = 0.5*(Sal(k-1) + Sal(k))
enddo
if (GV%Boussinesq .or. GV%semi_Boussinesq) then
call calculate_density_derivs(T_int, Sal_int, pressure, dbuoy_dT, dbuoy_dS, &
tv%eqn_of_state, (/2,nzc/), scale=-g_R0 )
else
! These should perhaps be combined into a single call to calculate the thermal expansion
! and haline contraction coefficients?
call calculate_specific_vol_derivs(T_int, Sal_int, pressure, dSpV_dT, dSpV_dS, &
tv%eqn_of_state, (/2,nzc/) )
call calculate_density(T_int, Sal_int, pressure, rho_int, tv%eqn_of_state, (/2,nzc/) )
do K=2,nzc
dbuoy_dT(K) = (US%L_to_Z**2 * GV%g_Earth) * (rho_int(K) * dSpV_dT(K))
dbuoy_dS(K) = (US%L_to_Z**2 * GV%g_Earth) * (rho_int(K) * dSpV_dS(K))
enddo
endif
elseif (GV%Boussinesq .or. GV%semi_Boussinesq) then
do K=1,nzc+1 ; dbuoy_dT(K) = -g_R0 ; dbuoy_dS(K) = 0.0 ; enddo
else
do K=1,nzc+1 ; dbuoy_dS(K) = 0.0 ; enddo
dbuoy_dT(1) = -(US%L_to_Z**2 * GV%g_Earth) / GV%Rlay(1)
do K=2,nzc
dbuoy_dT(K) = -(US%L_to_Z**2 * GV%g_Earth) / (0.5*(GV%Rlay(k-1) + GV%Rlay(k)))
enddo
dbuoy_dT(nzc+1) = -(US%L_to_Z**2 * GV%g_Earth) / GV%Rlay(nzc)
endif
! N2_debug(1) = 0.0 ; N2_debug(nzc+1) = 0.0
! do K=2,nzc
! N2_debug(K) = max((dbuoy_dT(K) * (T0xdz(k-1)*Idz(k-1) - T0xdz(k)*Idz(k)) + &
! dbuoy_dS(K) * (S0xdz(k-1)*Idz(k-1) - S0xdz(k)*Idz(k))) * &
! I_dz_int(K), 0.0)
! enddo
! This call just calculates N2 and S2.
call calculate_projected_state(kappa, u, v, T, Sal, 0.0, nzc, hlay, I_dz_int, dbuoy_dT, dbuoy_dS, &
CS%vel_underflow, u, v, T, Sal, N2, S2, GV, US)
! ----------------------------------------------------
! Iterate
! ----------------------------------------------------
dt_rem = dt
do K=1,nzc+1
K_Q(K) = 0.0
kappa_avg(K) = 0.0 ; tke_avg(K) = 0.0
local_src_avg(K) = 0.0
! Use the grid spacings to scale errors in the source.
if ( h_Int(K) > 0.0 ) &
local_src_avg(K) = 0.1 * k0dt * I_dz_int(K) / h_Int(K)
enddo
! call cpu_clock_end(id_clock_setup)
! do itt=1,CS%max_RiNo_it
do itt=1,CS%max_KS_it
! ----------------------------------------------------
! Calculate new values of u, v, rho, N^2 and S.
! ----------------------------------------------------
! call cpu_clock_begin(id_clock_KQ)
call find_kappa_tke(N2, S2, kappa, Idz, h_Int, dz_Int, dz_h_Int, I_L2_bdry, f2, &
nzc, CS, GV, US, K_Q, tke, kappa_out, kappa_src, local_src)
! call cpu_clock_end(id_clock_KQ)
! call cpu_clock_begin(id_clock_avg)
! Determine the range of non-zero values of kappa_out.
ks_kappa = GV%ke+1 ; ke_kappa = 0
do K=2,nzc ; if (kappa_out(K) > 0.0) then
ks_kappa = K ; exit
endif ; enddo
do k=nzc,ks_kappa,-1 ; if (kappa_out(K) > 0.0) then
ke_kappa = K ; exit
endif ; enddo
if (ke_kappa == nzc) kappa_out(nzc+1) = 0.0
! call cpu_clock_end(id_clock_avg)
! Determine how long to use this value of kappa (dt_now).
! call cpu_clock_begin(id_clock_project)
if ((ke_kappa < ks_kappa) .or. (itt==CS%max_KS_it)) then
dt_now = dt_rem
else
! Limit dt_now so that |k_src(k)-kappa_src(k)| < tol * local_src(k)
dt_test = dt_rem
do K=2,nzc
tol_max(K) = kappa_src(K) + tol_dksrc * local_src(K)
tol_min(K) = kappa_src(K) - tol_dksrc_low * local_src(K)
tol_chg(K) = tol2 * local_src_avg(K)
enddo
do itt_dt=1,(CS%max_KS_it+1-itt)/2
! The maximum number of times that the time-step is halved in
! seeking an acceptable timestep is reduced with each iteration,
! so that as the maximum number of iterations is approached, the
! whole remaining timestep is used. Typically, an acceptable
! timestep is found long before the minimum is reached, so the
! value of max_KS_it may be unimportant, especially if it is large
! enough.
call calculate_projected_state(kappa_out, u, v, T, Sal, 0.5*dt_test, nzc, hlay, I_dz_int, &
dbuoy_dT, dbuoy_dS, CS%vel_underflow, u_test, v_test, &
T_test, S_test, N2, S2, GV, US, ks_int=ks_kappa, ke_int=ke_kappa)
valid_dt = .true.
Idtt = 1.0 / dt_test
do K=max(ks_kappa-1,2),min(ke_kappa+1,nzc)
if (N2(K) < Ri_crit * S2(K)) then ! Equivalent to Ri < Ri_crit.
K_src(K) = (2.0 * CS%Shearmix_rate * sqrt(S2(K))) * &
((Ri_crit*S2(K) - N2(K)) / (Ri_crit*S2(K) + CS%FRi_curvature*N2(K)))
if (CS%restrictive_tolerance_check) then
if ((K_src(K) > min(tol_max(K), kappa_src(K) + Idtt*tol_chg(K))) .or. &