-
Notifications
You must be signed in to change notification settings - Fork 96
/
ocean_vert_kpp_mom4p1.F90
3314 lines (2823 loc) · 139 KB
/
ocean_vert_kpp_mom4p1.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
module ocean_vert_kpp_mom4p1_mod
!
!<CONTACT EMAIL="GFDL.Climate.Model.Info@noaa.gov"> A. Rosati
!</CONTACT>
!
!<CONTACT EMAIL="martin.schmidt@io-warnemuende.de"> Martin Schmidt
!</CONTACT>
!
!<REVIEWER EMAIL="wily@ucar.edu"> Bill Large
!</REVIEWER>
!
!<REVIEWER EMAIL="GFDL.Climate.Model.Info@noaa.gov"> Stephen Griffies
!</REVIEWER>
!
!<REVIEWER EMAIL="GFDL.Climate.Model.Info@noaa.gov"> M.J. Harrison
!</REVIEWER>
!
!<REVIEWER EMAIL="GFDL.Climate.Model.Info@noaa.gov"> Hyun-Chul Lee
!</REVIEWER>
!
!<OVERVIEW>
! Vertical viscosity and diffusivity according KPP.
! This module has extra code options to handle regions of extremely fresh water.
! This particular version of the KPP module is frozen based on its use in MOM4p1
! at GFDL for the IPCC AR5 climate model ESM2M.
!</OVERVIEW>
!
!<DESCRIPTION>
! This module computes vertical viscosity and diffusivity according to
! the K-profile parameterization scheme of Large, McWilliams, and
! Doney (1994). It computes both local and non-local mixing.
! The code has been updated to MOM4p1, so that vertical grid increments
! are suitable for generalized vertical coordinate models. When run
! as geopotential model, there will be some differences, since the
! MOM4.0 code (available in ocean_vert_kpp_mom4p0.F90) incorrectly
! ignored the free surface undulations affecting the top model grid
! cell thickness.
!
! This version of KPP has been implemented only for the Bgrid.
!</DESCRIPTION>
!
! <INFO>
!
! <REFERENCE>
! W.G. Large and J.C. McWilliams and S.C. Doney
! Oceanic vertical mixing: A review and a model with
! a nonlocal boundary layer parameterization
! Reviews of Geophysics (1994) vol 32 pages 363-403
! </REFERENCE>
! <REFERENCE>
! Danabasoglu etal (2006)
! Diurnal coupling in the tropical oceans of CCSM3
! Journal of Climate (2006) vol 19 pages 2347--2365
! </REFERENCE>
!
! <NOTE>
! Original numerical algorithm by Bill Large at NCAR June 6, 1994
! </NOTE>
!
! <NOTE>
! Equation numbers in the code refer to the Large etal paper.
! </NOTE>
!
! <NOTE>
! Surface fresh water contributes to surface buoyancy via conversion to
! a locally implied salt flux.
! </NOTE>
!
! </INFO>
!
!<NAMELIST NAME="ocean_vert_kpp_mom4p1_nml">
!
! <DATA NAME="use_this_module" TYPE="logical">
! Logical switch to enable kpp diffusion. Default is false.
! </DATA>
! <DATA NAME="debug_this_module" TYPE="logical">
! Logical switch for debugging. Default debug_this_module=.false.
! </DATA>
!
! <DATA NAME="shear_instability" TYPE="logical">
! logical switch for shear instability mixing.
! Default shear_instability=.true.
! </DATA>
! <DATA NAME="double_diffusion" TYPE="logical">
! Logical switch for double-diffusive mixing.
! Default double_diffusion=.true.
! </DATA>
!
! <DATA NAME="diff_cbt_iw" UNITS="m^2/sec" TYPE="real">
! Background vertical diffusivity. Note that if using Bryan-Lewis as a
! background diffusivity, then should set diff_cbt_iw=0.0.
! </DATA>
! <DATA NAME="visc_cbu_iw" UNITS="m^2/sec" TYPE="real">
! Background vertical viscosity
! </DATA>
! <DATA NAME="visc_cbu_limit" UNITS="m^2/sec" TYPE="real">
! Enhanced vertical viscosity due to shear instability
! </DATA>
! <DATA NAME="diff_cbt_limit" UNITS="m^2/sec" TYPE="real">
! Enhanced vertical diffusivity due to shear instability
! </DATA>
! <DATA NAME="visc_con_limit" UNITS="m^2/sec" TYPE="real">
! Enhanced vertical viscosity in regions of convection
! </DATA>
! <DATA NAME="diff_con_limit" UNITS="m^2/sec" TYPE="real">
! Enhanced vertical diffusivity in regions of convection
! </DATA>
!
! <DATA NAME="concv" UNITS="dimensionless" TYPE="real">
! constant for pure convection (eqn. 23 of Large etal)
! </DATA>
! <DATA NAME="Ricr" UNITS="dimensionless" TYPE="real">
! Critical bulk Richardson number. Default from NCAR is
! 0.3, though this number has a large uncertainty and some
! find that 1.0 can be of use.
! </DATA>
! <DATA NAME="non_local_kpp" TYPE="logical">
! logical switch for enabling the non-local mixing aspect of kpp.
! Default is .true. as this is what the original KPP scheme suggests.
! </DATA>
! <DATA NAME="smooth_blmc" TYPE="logical">
! Smooth boundary layer diffusitivies to remove grid scale noise.
! Such noise is apparent in the diagnosed mixed layer depth as well
! as the SST, especially when running coupled models where forcing
! has high temporal frequency.
! Default smooth_blmc=.false.
!
! Warning: This smoother can cause some problems with ghat in regions
! of zero surface forcing. To understand details, one needs
! the paper of Large et al. Vertical diffusion has the general form
! <wx> = K(x_z - ghats)
! In the surface layer a vertical scale function ws is estimated.
! We have K ~ ws and ghats ~1/ws. If wind stress is zero the vertical
! scale ws becomes zero too. Hence, ghats is very large
! (something finite, since it is divided by ws+epsln). Now it may happen,
! that the bouyancy flux becomes negative (~ -10-30). This enables
! the nonlocal scheme. Because the mixing coefficient in the
! surface boundary layer scales with ws the corresponding
! time tendency should be of the order (1/ws * ws = finite). However,
! if smooth_blmc is enabled, it may happen, that from neighbouring
! points with different mixing depth a finite value for
! the mixing coefficient leaks in. In this case
! the tracer time tendency from the nonlocal scheme becomes huge
! and the model fails.
!
! The smoother destroys the consistency between ghats and diff_cbt.
! In most cases this should not matter, but the example shows,
! that sudden model failure is possible under otherwise
! stable and smooth conditions.
!
! </DATA>
! <DATA NAME="kl_min" TYPE="integer">
! Lower loop index for finding new kbl. Needed for use with certain
! tests of OBC, where kl_min=1 needed, whereas default in original
! implementation has kl_min=2. Default in MOM is kl_min=2.
! </DATA>
! <DATA NAME="kbl_standard_method" TYPE="logical">
! For computing kbl as in the MOM4p0d code, which is taken from
! the original NCAR scheme. If false, then will slightly modify
! the logic. The modified logic has been found necessary when running
! with as few as two grid cells in the vertical.
! Default kbl_standard_method=.true.
! </DATA>
! <DATA NAME="limit_with_hekman" TYPE="logical">
! Limiting the boundary layer depth with the Ekman depth may result in a
! shallow boundary layer. In this case the internal values of the vertical
! mixing and viscosity coefficients may be large. This results in
! unrealistically large non-local vertical mixing
! Default limit_with_hekman=.true.
! </DATA>
! <DATA NAME="limit_ghats" TYPE="logical">
! Limits the non-local vertical tracer flux to the value of the tracer
! surface flux.
! Default limit_ghats=.false.
! </DATA>
! <DATA NAME="hbl_with_rit" TYPE="logical">
! The default method for determination of the boundary layer depth may fail
! if the water column is instable (negative Richardson number) below or above
! the layer that contains the diagnosed hbl.
! With hbl_with_rit=.true. the search for the boundary layer depth is continued
! downward in this case even if the bulk Richardson number exceeds the
! critical value. This removes a lot of noise from the boundary layer depth.
! Default hbl_with_rit=.false.
! </DATA>
! <DATA NAME="radiation_large" TYPE="logical">
! Remove the shortwave radiation leaving the boundary layer to the ocean interior
! (hence, not absorbed in the boundary layer) from non-local vertical heat flux
! Default radiation_large=.false.
! </DATA>
! <DATA NAME="radiation_zero" TYPE="logical">
! Remove the all shortwave radiation from non-local vertical heat flux.
! Default radiation_zero=.false.
! </DATA>
! <DATA NAME="radiation_iow" TYPE="logical">
! Keep only the shortwave radiation absorbed between the surface and a certain level
! in non-local vertical heat flux through this level.
! Default radiation_iow=.false.
! </DATA>
!
! <DATA NAME="bvf_from_below" TYPE="logical">
! Use BV-freq. at the cell bottom instead of the cell top
! as in Danabasoglu et al. (2006).
! Default bvf_from_below=.false., as this will recover
! older behaviour.
! </DATA>
! <DATA NAME="variable_vtc" TYPE="logical">
! Make vtc dependent on BV-freq. as in Danabasoglu et al. (2006).
! Default variable_vtc=.false., as this will recover
! older behaviour.
! </DATA>
! <DATA NAME="use_max_shear" TYPE="logical">
! Use maximum shear instead of 4-point average
! (as in Danabasoglu et al. (2006)).
! Default use_max_shear=.false., as this will recover
! older behaviour.
! </DATA>
!
! <DATA NAME="linear_hbl" TYPE="logical">
! Use linear interpolation to find the position of hbl.
! If set to false, then use the quadratic interpolation
! as in Danabasoglu et al. (2006). The quadratic approach
! generally yields a slightly deeper surface boundary layer.
! Default linear_hbl=.true., as this will recover
! older behaviour.
! </DATA>
! <DATA NAME="wsfc_combine_runoff_calve" TYPE="logical">
! For computing wsfc as in the MOM4p0d code, where we combine
! the runoff+calving into a single field called river.
! The alternative keeps the fields separate, as would be appropriate
! for a land model that separately tracks the tracer content in the
! calving and runoff.
! Default wsfc_combine_runoff_calve=.true., as this will recover
! the previous behaviour, to the bit.
! </DATA>
!
! <DATA NAME="smooth_ri_kmax_eq_kmu" TYPE="logical">
! When smoothing the Richardson number, we do so over a vertical
! column with max k-levels set by either kmt or kmu. The proper
! approach is kmu, since we are smoothing riu. But for backwards
! compatibility, we default to smooth_ri_kmax_eq_kmu=.false.
! </DATA>
!
!</NAMELIST>
use constants_mod, only: epsln, pi
use diag_manager_mod, only: register_diag_field, register_static_field
use fms_mod, only: FATAL, NOTE, WARNING, stdout, stdlog
use fms_mod, only: write_version_number, open_namelist_file, check_nml_error, close_file
use fms_mod, only: read_data
use mpp_domains_mod, only: mpp_update_domains, NUPDATE, EUPDATE, EDGEUPDATE
use mpp_mod, only: input_nml_file, mpp_error
use ocean_density_mod, only: density, density_delta_z, density_delta_sfc
use ocean_domains_mod, only: get_local_indices
use ocean_parameters_mod, only: grav, cp_ocean, missing_value, von_karman
use ocean_parameters_mod, only: ZSIGMA, PSIGMA, rho0r, rho0, PRESSURE_BASED
use ocean_types_mod, only: ocean_grid_type, ocean_domain_type
use ocean_types_mod, only: ocean_prog_tracer_type, ocean_diag_tracer_type
use ocean_types_mod, only: ocean_velocity_type, ocean_density_type
use ocean_types_mod, only: ocean_time_type, ocean_time_steps_type, ocean_thickness_type
use ocean_workspace_mod, only: wrk1, wrk2, wrk3, wrk4, wrk5
use ocean_workspace_mod, only: wrk1_2d, wrk2_2d
use ocean_util_mod, only: diagnose_2d, diagnose_3d, diagnose_sum
use ocean_tracer_util_mod, only: diagnose_3d_rho
implicit none
private
public ocean_vert_kpp_mom4p1_init
public vert_mix_kpp_mom4p1
private bldepth
private wscale
private ri_iwmix
private ddmix
private blmix_kpp
private enhance
private ri_for_kpp
private watermass_diag_init
private watermass_diag
private get_langmuir_number
private ust_2_u10_coare3p5
#include <ocean_memory.h>
#ifdef MOM_STATIC_ARRAYS
real, dimension(isd:ied,jsd:jed) :: bfsfc ! surface buoyancy forcing (m^2/s^3)
real, dimension(isd:ied,jsd:jed) :: ws ! scalar velocity scale (m/s)
real, dimension(isd:ied,jsd:jed) :: wm ! momentum velocity scale (m/s)
real, dimension(isd:ied,jsd:jed) :: Ustk2 ! magnitude of surface stokes drift velocity ^2 (m^2 / s^2)
real, dimension(isd:ied,jsd:jed) :: caseA ! = 1 in case A; =0 in case B
real, dimension(isd:ied,jsd:jed) :: stable ! = 1 in stable forcing; =0 in unstable
real, dimension(isd:ied,jsd:jed,3) :: dkm1 ! boundary layer diff_cbt at kbl-1 level
real, dimension(isd:ied,jsd:jed,nk,3) :: blmc ! boundary layer mixing coefficients
real, dimension(isd:ied,jsd:jed) :: sigma ! normalized depth (d / hbl)
real, dimension(isd:ied,jsd:jed) :: rhosfc ! potential density of sfc layer(kg/m^3)
real, dimension(isd:ied,jsd:jed,nk) :: talpha ! -d(rho)/ d(pot.temperature) (kg/m^3/C)
real, dimension(isd:ied,jsd:jed,nk) :: sbeta ! d(rho)/ d(salinity) (kg/m^3/PSU)
real, dimension(isd:ied,jsd:jed,nk) :: alphaDT ! alpha * DT across interfaces (kg/m^3)
real, dimension(isd:ied,jsd:jed,nk) :: betaDS ! beta * DS across interfaces (kg/m^3)
real, dimension(isd:ied,jsd:jed) :: Bo ! surface turb buoy. forcing (m^2/s^3)
real, dimension(isd:ied,jsd:jed) :: Bosol ! radiative buoy forcing (m^2/s^3)
real, dimension(isd:ied,jsd:jed,nk) :: dbloc ! local delta buoy at interfaces(m/s^2)
real, dimension(isd:ied,jsd:jed,nk) :: dbsfc ! delta buoy w/ respect to sfc (m/s^2)
real, dimension(isd:ied,jsd:jed,nk) :: dVsq ! (velocity shear re sfc)^2 (m/s)^2
real, dimension(isd:ied,jsd:jed,3) :: Rib ! Bulk Richardson number
real, dimension(isd:ied,jsd:jed,3) :: gat1
real, dimension(isd:ied,jsd:jed,3) :: dat1
type wsfc_type
real, dimension(isd:ied,jsd:jed) :: wsfc ! rho0r*(stf - pme*(t(i,j,k=1)-tpme) - river*(t(i,j,k=1)-triver))
end type wsfc_type
real, dimension(isd:ied,jsd:jed) :: sw_frac_hbl ! fractional shortwave penetration at base of mixed layer
integer, dimension(isd:ied,jsd:jed) :: kbl ! index of first grid level below hbl
real, private, dimension(isd:ied,jsd:jed,nk) :: riu ! Richardson number at base of U-cells
real, private, dimension(isd:ied,jsd:jed,nk) :: rit ! Richardson number at base of T-cells
real, private, dimension(isd:ied,jsd:jed,nk) :: ghats ! nonlocal transport (s/m^2)
real, private, dimension(isd:ied,jsd:jed) :: hblt ! boundary layer depth with tmask
real, private, dimension(isd:ied,jsd:jed) :: hbl ! boundary layer depth
real, dimension(isd:ied,jsd:jed) :: langmuir_factor ! Enhancnement due to langmuir turbulence3
real, dimension(isd:ied,jsd:jed) :: langmuir_number ! Langmuir number
#else
real, dimension(:,:), allocatable :: bfsfc ! surface buoyancy forcing (m^2/s^3)
real, dimension(:,:), allocatable :: ws ! scalar velocity scale (m/s)
real, dimension(:,:), allocatable :: wm ! momentum velocity scale (m/s)
real, dimension(:,:), allocatable :: Ustk2 ! Magnitude of surface stokes drift velocity ^2 (m^2/s^2)
real, dimension(:,:), allocatable :: caseA ! = 1 in case A; =0 in case B
real, dimension(:,:), allocatable :: stable ! = 1 in stable forcing; =0 in unstable
real, dimension(:,:,:), allocatable :: dkm1 ! boundary layer diff_cbt at kbl-1 level
real, dimension(:,:,:,:), allocatable :: blmc ! boundary layer mixing coefficients
real, dimension(:,:), allocatable :: sigma ! normalized depth (d / hbl)
real, dimension(:,:), allocatable :: rhosfc ! potential density of sfc layer(kg/m^3)
real, dimension(:,:,:), allocatable :: talpha ! -d(rho)/ d(pot.temperature) (kg/m^3/C)
real, dimension(:,:,:), allocatable :: sbeta ! d(rho)/ d(salinity) (kg/m^3/PSU)
real, dimension(:,:,:), allocatable :: alphaDT ! alpha * DT across interfaces (kg/m^3)
real, dimension(:,:,:), allocatable :: betaDS ! beta * DS across interfaces (kg/m^3)
real, dimension(:,:), allocatable :: Bo ! surface turb buoy. forcing (m^2/s^3)
real, dimension(:,:), allocatable :: Bosol ! radiative buoy forcing (m^2/s^3)
real, dimension(:,:,:), allocatable :: dbloc ! local delta buoy at interfaces(m/s^2)
real, dimension(:,:,:), allocatable :: dbsfc ! delta buoy w/ respect to sfc (m/s^2)
real, dimension(:,:,:), allocatable :: dVsq ! (velocity shear re sfc)^2 (m/s)^2
real, dimension(:,:,:), allocatable :: Rib ! Bulk Richardson number
real, dimension(:,:,:), allocatable :: gat1
real, dimension(:,:,:), allocatable :: dat1
type wsfc_type
real, dimension(:,:), pointer :: wsfc => NULL() ! rho0r*(stf - pme*(t(i,j,k=1)-tpme) - river*(t(i,j,k=1)-triver))
end type wsfc_type
real, dimension(:,:), allocatable :: sw_frac_hbl ! fractional shortwave penetration at base of mixed layer
integer, dimension(:,:), allocatable :: kbl ! index of first grid level below hbl
real, private, dimension(:,:,:), allocatable :: riu ! Richardson number at base of U-cells
real, private, dimension(:,:,:), allocatable :: rit ! Richardson number at base of T-cells
real, private, dimension(:,:,:), allocatable :: ghats ! nonlocal transport (s/m^2)
real, private, dimension(:,:), allocatable :: hblt ! boundary layer depth with tmask
real, private, dimension(:,:), allocatable :: hbl ! boundary layer depth
real, private, dimension(:,:), allocatable :: langmuir_factor ! Enhancnement due to langmuir turbulence3
real, private, dimension(:,:), allocatable :: langmuir_number ! Langmuir number
#endif
type(wsfc_type), dimension(:), allocatable :: wsfc
real, parameter :: epsilon = 0.1
real, parameter :: conc1 = 5.0
real, parameter :: zmin = -4.e-7 ! m3/s3 limit for lookup table of wm and ws
real, parameter :: zmax = 0.0 ! m3/s3 limit for lookup table of wm and ws
real, parameter :: umin = 0.0 ! m/s limit for lookup table of wm and ws
real, parameter :: umax = 0.04 ! m/s limit for lookup table of wm and ws
real :: Ricr = 0.3 ! critical bulk Richardson Number
real :: visc_cbu_limit = 50.0e-4 ! max visc due to shear instability
real :: diff_cbt_limit = 50.0e-4 ! max diff due to shear instability
real :: visc_con_limit = 0.1 ! m^2/s. visc due to convective instability
real :: diff_con_limit = 0.1 ! m^2/s. diff due to convective instability
real :: visc_cbu_iw = 1.0e-4 ! m^2/s. visc background due to internal waves
real :: diff_cbt_iw = 0.1e-4 ! m^2/s. diffusivity background due to internal waves
real :: Vtc ! non-dimensional coefficient for velocity
! scale of turbulant velocity shear
! (=function of concv,concs,epsilon,von_karman,Ricr)
real :: cg ! non-dimensional coefficient for counter-gradient term
real :: deltaz ! delta zehat in table
real :: deltau ! delta ustar in table
real :: concv = 1.8 ! constant for pure convection (eqn. 23)
real :: concv_r ! inverse concv
real :: vtc_flag = 0.0 ! default to the older approach.
real :: Lgam = 1.04 ! adjustment to non-gradient flux (McWilliam & Sullivan 2000)
real :: Cw_0 = 0.15 ! eq. (13) in Smyth et al (2002)
real :: l_smyth = 2.0 ! eq. (13) in Smyth et al (2002)
real :: LTmax = 5.0 ! maximum Langmuir turbulence enhancement factor (langmuir_factor) allowed
real :: Wstfac = 0.6 ! stability adjustment coefficient, eq. (13) in Smyth et al (2002)
! for global area normalization
real :: cellarea_r
! for vertical coordinate
integer :: vert_coordinate
integer :: vert_coordinate_class
real :: rho_cp
real :: inv_rho_cp
integer :: kl_min=2
logical :: kbl_standard_method=.true.
logical :: use_this_module = .false. ! logical switch for turning on the kpp scheme
logical :: shear_instability = .true. ! logical switch for shear instability mixing
logical :: double_diffusion = .true. ! logical switch for double-diffusive mixing
logical :: limit_ghats = .false. ! logical switch to limit ghats*diff to 1
logical :: limit_with_hekman = .true. ! logical switch to limit hbl with hekman
logical :: radiation_large = .false. ! logical switch to enable short wave radiation in non-local scheme h_gamma=h
logical :: radiation_zero = .false. ! logical switch to enable short wave radiation in non-local scheme h_gamma=0
logical :: radiation_iow = .false. ! logical switch to enable short wave radiation in non-local scheme IOW-method
logical :: hbl_with_rit = .false. ! logical switch to enable hbl correction for negative rit
logical :: use_sbl_bottom_flux = .false. ! logical switch to add flux through the sbl bottom to the non-local term
logical :: wsfc_combine_runoff_calve = .true. ! for combining runoff+calving to compute wfsc
logical :: bvf_from_below = .false. ! Use BV-freq. at the cell bottom instead of the cell top (Danabasoglu et al.)
logical :: variable_vtc = .false. ! Make vtc dependent on BV-freq.
logical :: use_max_shear = .false. ! Use maximum shear instead of 4-point average
logical :: linear_hbl = .true. ! To use the linear interpolation as Large etal (1994).
! Set to .false. for quadratic interpolation as in Danabasoglu et al.
logical :: calc_visc_on_cgrid=.false. ! calculate viscosity directly on c-grid
logical :: smooth_ri_kmax_eq_kmu=.false. ! to set details for smoothing the richardson number
real :: shear_instability_flag = 1.0 ! set to 1.0 if shear_instability=.true.
! internally set for computing watermass diagnstics
logical :: compute_watermass_diag = .false.
! for diagnostics
integer, dimension(:), allocatable :: id_nonlocal(:)
integer, dimension(:), allocatable :: id_nonlocal_on_nrho(:)
integer, dimension(:), allocatable :: id_ghats(:)
integer, dimension(:), allocatable :: id_wsfc(:)
integer, dimension(:), allocatable :: id_wbot(:)
logical :: used
integer :: id_diff_cbt_kpp_t =-1
integer :: id_diff_cbt_kpp_s =-1
integer :: id_diff_cbt_conv =-1
integer :: id_hblt =-1
integer :: id_ws =-1
integer :: id_lang_enh =-1
integer :: id_lang =-1
integer :: id_u10 =-1
integer :: id_neut_rho_kpp_nloc =-1
integer :: id_pot_rho_kpp_nloc =-1
integer :: id_wdian_rho_kpp_nloc =-1
integer :: id_tform_rho_kpp_nloc =-1
integer :: id_neut_rho_kpp_nloc_on_nrho =-1
integer :: id_wdian_rho_kpp_nloc_on_nrho =-1
integer :: id_tform_rho_kpp_nloc_on_nrho =-1
integer :: id_eta_tend_kpp_nloc =-1
integer :: id_eta_tend_kpp_nloc_glob=-1
integer :: id_neut_temp_kpp_nloc =-1
integer :: id_wdian_temp_kpp_nloc =-1
integer :: id_tform_temp_kpp_nloc =-1
integer :: id_neut_temp_kpp_nloc_on_nrho =-1
integer :: id_wdian_temp_kpp_nloc_on_nrho =-1
integer :: id_tform_temp_kpp_nloc_on_nrho =-1
integer :: id_neut_salt_kpp_nloc =-1
integer :: id_wdian_salt_kpp_nloc =-1
integer :: id_tform_salt_kpp_nloc =-1
integer :: id_neut_salt_kpp_nloc_on_nrho =-1
integer :: id_wdian_salt_kpp_nloc_on_nrho =-1
integer :: id_tform_salt_kpp_nloc_on_nrho =-1
logical :: non_local_kpp = .true. ! enable/disable non-local term in KPP
logical :: smooth_blmc = .false. ! smooth boundary layer diffusitivies to remove grid scale noise
logical :: do_langmuir = .false. ! whether or not calcualte langmuir turbulence enhance factor
logical :: do_langmuir_cvmix = .false. ! whether or not calcualte langmuir turbulence enhance factor via cvmix method based on Van Roekel 2012
logical :: calculate_u10 = .false. ! If 10m not available the estimate u10 from u_star
integer, parameter :: nni = 890 ! number of values for zehat in the look up table
integer, parameter :: nnj = 480 ! number of values for ustar in the look up table
real, dimension(0:nni+1,0:nnj+1) :: wmt ! lookup table for wm, the turbulent velocity scale for momentum
real, dimension(0:nni+1,0:nnj+1) :: wst ! lookup table for ws, the turbulent velocity scale scalars
real :: tracer_timestep = 0
type(ocean_domain_type), pointer :: Dom => NULL()
type(ocean_grid_type), pointer :: Grd => NULL()
integer :: num_prog_tracers=0, index_temp, index_salt
integer :: num_diag_tracers=0, index_frazil
integer :: kbl_max=2 ! helps to limit vertikal loops
character(len=256) :: version=&
'$Id: ocean_vert_kpp_mom4p1.F90,v 20.0 2013/12/14 00:16:44 fms Exp $'
character (len=128) :: tagname = &
'$Name: tikal $'
logical :: module_is_initialized = .FALSE.
logical :: debug_this_module = .FALSE.
namelist /ocean_vert_kpp_mom4p1_nml/ use_this_module, shear_instability, double_diffusion, &
diff_cbt_iw, visc_cbu_iw, &
visc_cbu_limit, diff_cbt_limit, &
visc_con_limit, diff_con_limit, &
concv, Ricr, non_local_kpp, smooth_blmc, &
Lgam, Cw_0,l_smyth, LTmax, Wstfac, &
kl_min, kbl_standard_method, debug_this_module, &
limit_with_hekman, limit_ghats, hbl_with_rit, &
radiation_large, radiation_zero, radiation_iow, &
use_sbl_bottom_flux, wsfc_combine_runoff_calve, &
bvf_from_below, variable_vtc, use_max_shear, &
linear_hbl, calc_visc_on_cgrid, smooth_ri_kmax_eq_kmu, &
do_langmuir, do_langmuir_cvmix, calculate_u10
contains
!#######################################################################
! <SUBROUTINE NAME="ocean_vert_kpp_mom4p1_init">
!
! <DESCRIPTION>
! Initialization for the KPP vertical mixing scheme
! </DESCRIPTION>
!
subroutine ocean_vert_kpp_mom4p1_init (Grid, Domain, Time, Time_steps, Dens, T_prog, T_diag, &
ver_coordinate, ver_coordinate_class)
type(ocean_grid_type), intent(in), target :: Grid
type(ocean_domain_type), intent(in), target :: Domain
type(ocean_time_type), intent(in) :: Time
type(ocean_time_steps_type), intent(in) :: Time_steps
type(ocean_density_type), intent(in) :: Dens
type(ocean_prog_tracer_type), intent(inout) :: T_prog(:)
type(ocean_diag_tracer_type), intent(in) :: T_diag(:)
integer, intent(in) :: ver_coordinate
integer, intent(in) :: ver_coordinate_class
real, parameter :: cstar = 10.0 ! proportionality coefficient for nonlocal transport
real, parameter :: conam = 1.257
real, parameter :: concm = 8.380
real, parameter :: conc2 = 16.0
real, parameter :: zetam = -0.2
real, parameter :: conas = -28.86
real, parameter :: concs = 98.96
real, parameter :: conc3 = 16.0
real, parameter :: zetas = -1.0
real :: zehat ! = zeta * ustar**3
real :: zeta ! = stability parameter d/L
real :: usta
integer :: i, j, n, ioun, io_status, ierr
integer :: rib_dim = 3
integer :: stdoutunit,stdlogunit
stdoutunit=stdout();stdlogunit=stdlog()
if ( module_is_initialized ) then
call mpp_error(FATAL, &
'==>Error from ocean_vert_kpp_mom4p1_mod (ocean_vert_kpp_mom4p1_init): module is already initialized')
endif
module_is_initialized = .TRUE.
vert_coordinate = ver_coordinate
vert_coordinate_class = ver_coordinate_class
rho_cp = rho0*cp_ocean
inv_rho_cp = 1.0/rho_cp
call write_version_number(version, tagname)
! provide for namelist over-ride of defaults
#ifdef INTERNAL_FILE_NML
read (input_nml_file, nml=ocean_vert_kpp_mom4p1_nml, iostat=io_status)
ierr = check_nml_error(io_status,'ocean_vert_kpp_mom4p1_nml')
#else
ioun = open_namelist_file ()
read (ioun, ocean_vert_kpp_mom4p1_nml,iostat=io_status)
ierr = check_nml_error(io_status,'ocean_vert_kpp_mom4p1_nml')
call close_file(ioun)
#endif
write (stdoutunit,'(/)')
write (stdoutunit, ocean_vert_kpp_mom4p1_nml)
write (stdlogunit, ocean_vert_kpp_mom4p1_nml)
#ifndef MOM_STATIC_ARRAYS
call get_local_indices(Domain,isd,ied,jsd,jed,isc,iec,jsc,jec)
nk = Grid%nk
#endif
kbl_max=nk
if(use_this_module) then
write(stdoutunit,'(/1x,a)') '==> NOTE: USING KPP vertical mixing scheme.'
write(stdoutunit,'(1x,a)') '==> NOTE: KPP is typically run with penetrative shortwave heating.'
write(stdoutunit,'(1x,a/)') '==> NOTE: KPP is typically run with a seasonal and/or diurnal cycle.'
else
write(stdoutunit,'(/1x,a)')'==> NOTE: NOT USING KPP vertical mixing scheme.'
return
endif
if(debug_this_module) then
write(stdoutunit,'(/1x,a)')'==> NOTE: Running KPP with debug_this_module=.true. Set vmixing coeffs to constant.'
endif
if(vert_coordinate==ZSIGMA .or. vert_coordinate==PSIGMA) then
call mpp_error(WARNING,&
'==>ocean_vert_kpp_mom4p1_mod: ZSIGMA or PSIGMA are not fully working w/ KPP. ')
call mpp_error(WARNING,&
'Problems exist /w tracer balances not closing AND instabilities w/ nonlocal transport.')
endif
write(stdoutunit,'(/a,f10.2)') &
'==>Note from ocean_vert_kpp_mom4p1_mod: using forward time step for vert-frict of (secs)', &
Time_steps%dtime_u
write(stdoutunit,'(/a,f10.2)') &
'==>Note from ocean_vert_kpp_mom4p1_mod: using forward time step for vert-diff of (secs)', &
Time_steps%dtime_t
if(kbl_standard_method) then
write(stdoutunit,'(1x,a)') '==> NOTE from ocean_vert_kpp_mom4p1_mod: adjust kbl to hbl with the standard method.'
else
write(stdoutunit,'(1x,a)') '==> NOTE from ocean_vert_kpp_mom4p1_mod: adjust kbl to hbl with the non-standard method.'
if (kl_min > 1) write(stdoutunit,'(1x,a,I4)') '==>WARNING from ocean_vert_kpp_mom4p1_mod:'// &
' change kl_min to 1 to avoid negative mixing coefficients with low wind stress.'
endif
if(use_sbl_bottom_flux) then
write(stdoutunit,'(1x,a)') '==> NOTE from ocean_vert_kpp_mom4p1_mod: '// &
'Add tracer flux through the SBL-bottom to the nonlocal scheme. This is experimental.'
endif
if(bvf_from_below) then
write(stdoutunit,'(1x,a)') '==> NOTE from ocean_vert_kpp_mom4p1_mod: '// &
'Find BVF from a forward derivatives.'
else
write(stdoutunit,'(1x,a)') '==> NOTE from ocean_vert_kpp_mom4p1_mod: '// &
'Find BVF from a backward derivatives.'
endif
if(variable_vtc) then
write(stdoutunit,'(1x,a)') '==> NOTE from ocean_vert_kpp_mom4p1_mod: '// &
'Diagnose hbl with concv that depends on BVF.'
else
write(stdoutunit,'(1x,a)') '==> NOTE from ocean_vert_kpp_mom4p1_mod: '// &
'Diagnose hbl with constant concv.'
endif
if(use_max_shear) then
write(stdoutunit,'(1x,a)') '==> NOTE from ocean_vert_kpp_mom4p1_mod: '// &
'Use maximum shear around a tracer cell for diagnostics of hbl.'
else
write(stdoutunit,'(1x,a)') '==> NOTE from ocean_vert_kpp_mom4p1_mod: '// &
'Use average shear around a tracer cell for diagnostics of hbl.'
endif
if (linear_hbl) then
write(stdoutunit,'(1x,a)') '==> NOTE from ocean_vert_kpp_mom4p1_mod: Use linear interpolation to find hbl'
else
write(stdoutunit,'(1x,a)') '==> NOTE from ocean_vert_kpp_mom4p1_mod: Use quadratic interpolation to find hbl'
endif
if(limit_ghats) then
write(stdoutunit,'(1x,a)') '==> NOTE from ocean_vert_kpp_mom4p1_mod: Limit ghats*diff_cbt to 1.'
else
write(stdoutunit,'(1x,a)') '==> NOTE from ocean_vert_kpp_mom4p1_mod: Do not limit ghats*diff_cbt to 1.'
endif
if(limit_with_hekman) then
write(stdoutunit,'(1x,a)') '==> NOTE from ocean_vert_kpp_mom4p1_mod: Limit hbl with hekman for stable case.'
else
write(stdoutunit,'(1x,a)') '==> NOTE from ocean_vert_kpp_mom4p1_mod: Do not limit hbl with hekman for stable case.'
endif
if(calc_visc_on_cgrid) then
write(stdoutunit,'(1x,a)') '==> NOTE from ocean_vert_kpp_mom4p1_mod: Calculate viscosity at c-grid.'
else
write(stdoutunit,'(1x,a)') '==> NOTE from ocean_vert_kpp_mom4p1_mod: Interpolate viscosity to c-grid.'
endif
if(radiation_large) then
write(stdoutunit,'(1x,a)') '==> NOTE from ocean_vert_kpp_mom4p1_mod: sw-radiation in non-local scheme h_gamma=h.'
elseif (radiation_zero) then
write(stdoutunit,'(1x,a)') '==> NOTE from ocean_vert_kpp_mom4p1_mod: sw-radiation in non-local scheme with h_gamma=0.'
elseif (radiation_iow) then
write(stdoutunit,'(1x,a)') '==> NOTE from ocean_vert_kpp_mom4p1_mod: sw-radiation in non-local scheme with IOW-method.'
else
write(stdoutunit,'(1x,a)') '==> NOTE from ocean_vert_kpp_mom4p1_mod: Leave full sw-radiation in non-local surface flux.'
endif
if(do_langmuir .and. .not. do_langmuir_cvmix) then
write(stdoutunit,'(1x,a)') '==> NOTE from ocean_vert_kpp_mom4p1_mod: Enhancing mixing via langmuir turbulence according to Smyth 2002.'
endif
if(do_langmuir .and. do_langmuir_cvmix) then
write(stdoutunit,'(1x,a)') '==> NOTE from ocean_vert_kpp_mom4p1_mod: Enhancing mixing via langmuir turbulence according to van Roekel 2012.'
endif
if( calculate_u10 ) then
write(stdoutunit,'(1x,a)') '==> NOTE from ocean_vert_kpp_mom4p1_mod: Calculating U_10 on T-grid from ustar COARE 3.5 paper (Edson et al., 2013).'
endif
if(radiation_large .and. radiation_zero) call mpp_error(FATAL,&
'==>ocean_vert_kpp_mom4p1_mod: Do not enable radiation_large and radiation_zero together. ')
if(radiation_large .and. radiation_iow) call mpp_error(FATAL,&
'==>ocean_vert_kpp_mom4p1_mod: Do not enable radiation_large and radiation_iow together. ')
if(radiation_zero .and. radiation_iow) call mpp_error(FATAL,&
'==>ocean_vert_kpp_mom4p1_mod: Do not enable radiation_zero and radiation_iow together. ')
tracer_timestep = Time_steps%dtts
index_temp=-1;index_salt=-1
num_prog_tracers = size(T_prog(:))
do n= 1, num_prog_tracers
if (T_prog(n)%name == 'temp') index_temp = n
if (T_prog(n)%name == 'salt') index_salt = n
enddo
allocate( wsfc(num_prog_tracers) )
index_frazil=-1
num_diag_tracers = size(T_diag(:))
do n= 1, num_diag_tracers
if (T_diag(n)%name == 'frazil') index_frazil = n
enddo
if (index_temp < 1 .or. index_salt < 1) then
call mpp_error(FATAL, &
'==>Error: ocean_vert_kpp_mom4p1_mod: temp and/or salt not present in tracer array')
endif
! for computing Vtsq
if (variable_vtc) then
vtc_flag=1.0 ! as in Danabasoglu etal (2006)
else
vtc_flag=0.0 ! as in Large etal (1994)
endif
Dom => Domain
Grd => Grid
if (linear_hbl) rib_dim=2
cellarea_r = 1.0/(epsln + Grd%tcellsurf)
#ifndef MOM_STATIC_ARRAYS
if(radiation_iow) then
allocate( T_prog(index_temp)%radiation(isd:ied,jsd:jed,nk))
T_prog(index_temp)%radiation(:,:,:) = 0.0
endif
allocate (ghats(isd:ied,jsd:jed,nk))
allocate (riu(isd:ied,jsd:jed,nk))
allocate (rit(isd:ied,jsd:jed,nk))
allocate (hblt(isd:ied,jsd:jed))
allocate (hbl(isd:ied,jsd:jed))
do n=1,num_prog_tracers
allocate( wsfc(n)%wsfc(isd:ied,jsd:jed) )
enddo
allocate (bfsfc(isd:ied,jsd:jed)) ! surface buoyancy forcing (m^2/s^3)
allocate (caseA(isd:ied,jsd:jed)) ! = 1 in case A; =0 in case B
allocate (stable(isd:ied,jsd:jed)) ! = 1 in stable forcing; =0 in unstable
allocate (dkm1(isd:ied,jsd:jed,3)) ! boundary layer diff_cbt at kbl-1 level
allocate (blmc(isd:ied,jsd:jed,nk,3)) ! boundary layer mixing coefficients
allocate (sigma(isd:ied,jsd:jed)) ! normalized depth (d / hbl)
allocate (rhosfc(isd:ied,jsd:jed)) ! potential density of sfc layer(g/m^3)
allocate (talpha(isd:ied,jsd:jed,nk)) ! -d(rho)/ d(pot.temperature) (g/m^3/C)
allocate (sbeta(isd:ied,jsd:jed,nk)) ! d(rho)/ d(salinity) (g/m^3/PSU)
allocate (alphaDT(isd:ied,jsd:jed,nk)) ! alpha * DT across interfaces (g/m^3)
allocate (betaDS(isd:ied,jsd:jed,nk)) ! beta * DS across interfaces (g/m^3)
allocate (Bo(isd:ied,jsd:jed)) ! surface turb buoy. forcing (m^2/s^3)
allocate (Bosol(isd:ied,jsd:jed)) ! radiative buoy forcing (m^2/s^3)
allocate (dbloc(isd:ied,jsd:jed,nk)) ! local delta buoy at interfaces(m/s^2)
allocate (dbsfc(isd:ied,jsd:jed,nk)) ! delta buoy w/ respect to sfc (m/s^2)
allocate (dVsq(isd:ied,jsd:jed,nk)) ! (velocity shear re sfc)^2 (m/s)^2
allocate (kbl(isd:ied,jsd:jed)) ! index of first grid level below hbl
allocate (Rib(isd:ied,jsd:jed,rib_dim)) ! Bulk Richardson number
allocate (wm(isd:ied,jsd:jed)) ! momentum turbulent velocity scales (m/s)
allocate (ws(isd:ied,jsd:jed)) ! scalar turbulent velocity scales (m/s)
allocate (Ustk2(isd:ied,jsd:jed)) ! Magnitude of surface stokes drift velocity ^2 (m^2/s^2)
allocate (gat1(isd:ied,jsd:jed,3))
allocate (dat1(isd:ied,jsd:jed,3))
allocate(sw_frac_hbl(isd:ied,jsd:jed))
allocate (langmuir_factor(isd:ied,jsd:jed))
allocate (langmuir_number(isd:ied,jsd:jed))
#endif
kbl(:,:) = 0
ghats(:,:,:) = 0.0
riu(:,:,:) = 0.0
rit(:,:,:) = 0.0
hblt(:,:) = 0.0
hbl(:,:) = 0.0
sw_frac_hbl(:,:) = 0.0
Ustk2(:,:) = 0.0
do n = 1, num_prog_tracers
wsfc(n)%wsfc(:,:) = 0.0
enddo
!-----------------------------------------------------------------------
! initialize some constants for kmix subroutines, and initialize
! for kmix subroutine "wscale" the 2D-lookup table for wm and ws
! as functions of ustar and zetahat (=von_karman*sigma*hbl*bfsfc).
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
! define some non-dimensional constants (recall epsilon=0.1)
!-----------------------------------------------------------------------
concv_r = 1.0/concv
! Vtc used in eqn. 23
Vtc = concv * sqrt(0.2/concs/epsilon) / von_karman**2 / Ricr
! cg = cs in eqn. 20
cg = cstar * von_karman * (concs * von_karman * epsilon)**(1./3.)
!-----------------------------------------------------------------------
! construct the wm and ws lookup tables (eqn. 13 & B1)
!-----------------------------------------------------------------------
deltaz = (zmax-zmin)/(nni+1)
deltau = (umax-umin)/(nnj+1)
do i=0,nni+1
zehat = deltaz*(i) + zmin
do j=0,nnj+1
usta = deltau*(j) + umin
zeta = zehat/(usta**3+epsln)
if(zehat >= 0.) then
wmt(i,j) = von_karman*usta/(1.+conc1*zeta)
wst(i,j) = wmt(i,j)
else
if(zeta > zetam) then
wmt(i,j) = von_karman* usta * (1.-conc2*zeta)**(1./4.)
else
wmt(i,j) = von_karman* (conam*usta**3-concm*zehat)**(1./3.)
endif
if(zeta > zetas) then
wst(i,j) = von_karman* usta * (1.-conc3*zeta)**(1./2.)
else
wst(i,j) = von_karman* (conas*usta**3-concs*zehat)**(1./3.)
endif
endif
enddo
enddo
! set flag for shear instability mixing
if (shear_instability) then
write (stdoutunit,'(a)') 'Computing vertical mixing from shear instability in KPP module.'
shear_instability_flag=1.0
else
write (stdoutunit,'(a)') 'NOT computing vertical mixing from shear instability in KPP module.'
shear_instability_flag=0.0
endif
!-----------------------------------------------------------------------
! error checks and diagnostics setup
!-----------------------------------------------------------------------
if (Time_steps%aidif /= 1.0) then
call mpp_error(FATAL,&
'==>Error in ocean_vert_kpp_mom4p1_mod (ocean_vert_kpp_mom4p1_init): KPP must use aidif=1 for implicit vertical mixing')
endif
! register diagnostics
allocate(id_nonlocal(num_prog_tracers))
allocate(id_nonlocal_on_nrho(num_prog_tracers))
allocate(id_ghats(2))
allocate(id_wsfc(num_prog_tracers))
allocate(id_wbot(num_prog_tracers))
id_nonlocal=-1
id_nonlocal_on_nrho=-1
do n = 1, num_prog_tracers
if(n==index_temp) then
id_nonlocal(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_nonlocal_KPP', &
Grd%tracer_axes(1:3), Time%model_time, &
'cp*rho*dzt*nonlocal tendency from KPP', trim(T_prog(n)%flux_units), &
missing_value=missing_value, range=(/-1.e10,1.e10/))
id_nonlocal_on_nrho(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_nonlocal_KPP_on_nrho', &
Dens%neutralrho_axes(1:3), Time%model_time, &
'cp*rho*dzt*nonlocal tendency from KPP binned to neutral density', trim(T_prog(n)%flux_units), &
missing_value=missing_value, range=(/-1.e20,1.e20/))
id_wsfc(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_wsfc_KPP', &
Grd%tracer_axes(1:2), Time%model_time, &
'cp*rho*dzt*surface tendency from KPP', trim(T_prog(n)%flux_units), &
missing_value=missing_value, range=(/-1.e10,1.e10/))
else
id_nonlocal(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_nonlocal_KPP', &
Grd%tracer_axes(1:3), Time%model_time, &
'rho*dzt*nonlocal tendency from KPP', trim(T_prog(n)%flux_units), &
missing_value=missing_value, range=(/-1.e10,1.e10/))
id_nonlocal_on_nrho(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_nonlocal_KPP_on_nrho', &
Dens%neutralrho_axes(1:3), Time%model_time, &
'rho*dzt*nonlocal tendency from KPP binned to neutral density', trim(T_prog(n)%flux_units), &
missing_value=missing_value, range=(/-1.e20,1.e20/))
id_wsfc(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_wsfc_KPP', &
Grd%tracer_axes(1:2), Time%model_time, &
'rho*dzt*surface tendency from KPP', trim(T_prog(n)%flux_units), &
missing_value=missing_value, range=(/-1.e10,1.e10/))
endif
id_wbot(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_wbot_KPP', &
Grd%tracer_axes(1:2), Time%model_time, &
'tracer flux through sbl-bottom', trim(T_prog(n)%flux_units), &
missing_value=missing_value, range=(/-1.e10,1.e10/))
enddo
id_ghats(1) = register_diag_field ('ocean_model', 'temp_ghats_KPP', &
Grd%tracer_axes(1:3), Time%model_time, &
'nonlocal term ghats * diff_cbt from KPP', 'none', &
missing_value=missing_value, range=(/-1.e10,1.e10/))
id_ghats(2) = register_diag_field ('ocean_model', 'salt_ghats_KPP', &
Grd%tracer_axes(1:3), Time%model_time, &
'nonlocal term ghats * diff_cbt from KPP', 'none', &
missing_value=missing_value, range=(/-1.e10,1.e10/))
id_diff_cbt_kpp_t = register_diag_field('ocean_model','diff_cbt_kpp_t', &
Grd%tracer_axes(1:3),Time%model_time, 'vert diffusivity from kpp for temp',&
'm^2/sec', missing_value = missing_value, range=(/-1.e5,1.e5/))
id_diff_cbt_kpp_s = register_diag_field('ocean_model','diff_cbt_kpp_s', &
Grd%tracer_axes(1:3), Time%model_time, 'vert diffusivity from kpp for salt',&
'm^2/sec', missing_value = missing_value, range=(/-1.e5,1.e5/))
id_diff_cbt_conv = register_diag_field('ocean_model','diff_cbt_conv', &
Grd%tracer_axes(1:3),Time%model_time, 'vert diffusivity from kpp convection',&
'm^2/sec', missing_value = missing_value, range=(/-1.e5,1.e5/))
id_hblt = register_diag_field('ocean_model','hblt',Grd%tracer_axes(1:2), &
Time%model_time, 'T-cell boundary layer depth from KPP', 'm', &
missing_value = missing_value, range=(/-1.e5,1.e6/), &
standard_name='ocean_mixed_layer_thickness_defined_by_mixing_scheme')
id_ws = register_diag_field('ocean_model','wscale',Grd%tracer_axes(1:3), &
Time%model_time, 'wscale from KPP', 'm', &
missing_value = missing_value, range=(/-1.e5,1.e6/))
id_lang_enh = register_diag_field('ocean_model','lang_enh',Grd%tracer_axes(1:2), &
Time%model_time, 'enhancement due to langmuir turbulence', 'none', &
missing_value = missing_value, range=(/0.0,10./))
id_lang = register_diag_field('ocean_model','lang_num',Grd%tracer_axes(1:2), &
Time%model_time, 'Langmuir number', 'none', &
missing_value = missing_value, range=(/0.0,1.E10/))
id_u10 = register_diag_field('ocean_model','u10',Grd%tracer_axes(1:2), &
Time%model_time, '10m wind speed used for kpp Langmuir turbulence', 'm/s', &
missing_value = missing_value, range=(/0.0,1.e3/))
call watermass_diag_init(Time,Dens)
end subroutine ocean_vert_kpp_mom4p1_init
! </SUBROUTINE> NAME="ocean_vert_kpp_mom4p1_init">
!#######################################################################
! <SUBROUTINE NAME="vert_mix_kpp_mom4p1">
!
! <DESCRIPTION>
! This subroutine computes the vertical diffusivity and viscosity according
! to the KPP scheme of Large etal. In brief, the scheme does the
! following:
!
! --Compute interior mixing everywhere:
! interior mixing gets computed at all cell interfaces due to constant
! internal wave background activity ("visc_cbu_iw" and "diff_cbt_iw").
! Mixing is enhanced in places of static instability (local Ri < 0).
! Additionally, mixing can be enhanced by contribution from shear
! instability which is a function of the local Ri.
!
! --Double diffusion:
! Interior mixing can be enhanced by double diffusion due to salt
! fingering and diffusive convection ("double_diffusion=.true.").
!
! --Boundary layer:
!
! (A) Boundary layer depth:
! at every gridpoint the depth of the oceanic boundary layer
! ("hbl") gets computed by evaluating bulk richardson numbers.
!
! (B) Boundary layer mixing:
! within the boundary layer, above hbl, vertical mixing is
! determined by turbulent surface fluxes, and interior mixing at
! the lower boundary, i.e. at hbl.
!
! </DESCRIPTION>
!
subroutine vert_mix_kpp_mom4p1 (aidif, Time, Thickness, Velocity, T_prog, T_diag, Dens, &
swflx, sw_frac_zt, pme, river, visc_cbu, diff_cbt, &
diff_cbt_conv, hblt_depth, do_wave)
real, intent(in) :: aidif