-
Notifications
You must be signed in to change notification settings - Fork 96
/
ocean_submesoscale.F90
4469 lines (3762 loc) · 174 KB
/
ocean_submesoscale.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_submesoscale_mod
#define COMP isc:iec,jsc:jec
!
!<CONTACT EMAIL="[email protected]"> S. M. Griffies
!</CONTACT>
!
!<OVERVIEW>
! This module computes a streamfunction within
! the upper surface boundary layer, and applies this
! streamfunction to all tracers. It also optionally
! applies horizontal diffusion in the surface layer
! as determined by the strength of the streamfunction.
!</OVERVIEW>
!
!<DESCRIPTION>
! This module computes a streamfunction within
! the upper surface boundary layer, and applies this
! streamfunction to all tracers. It also optionally
! applies horizontal diffusion in the surface layer
! as determined by the strength of the streamfunction.
!</DESCRIPTION>
!
! <INFO>
!
! <REFERENCE>
! Fox-Kemper, Ferrari, and Hallberg 2008: Parameterization of
! mixed layer eddies. Part I: theory and diagnosis
! Journal of Physical Oceanography, vol. 38, pages 1145-1165.
! </REFERENCE>
!
! <REFERENCE>
! Fox-Kemper, Danabasoglu, Ferrari, and Hallberg 2008:
! Parameterizing submesoscale physics in global models.
! Clivar Exchanges, vol 13, no.1, Jan2008. pages 3-5.
! </REFERENCE>
!
! <REFERENCE>
! Fox-Kemper, Danabasoglu, Ferrari, Griffies, Hallberg,
! Holland, Peacock, Samuels, 2011: Parameterization of
! Mixed Layer Eddies. III: Global Implementation and
! Impact on Ocean Climate Simulations, Ocean Modelling,
! vol. 39, pages 61-78.
! </REFERENCE>
!
! <REFERENCE>
! Griffies, 2012: Elements of MOM
! </REFERENCE>
!
! </INFO>
!
!<NAMELIST NAME="ocean_submesoscale_nml">
! <DATA NAME="use_this_module=" TYPE="logical">
! Must be .true. to use this module.
! </DATA>
! <DATA NAME="debug_this_module" TYPE="logical">
! For debugging purposes.
! </DATA>
! <DATA NAME="diag_step" TYPE="integer">
! Number of time steps between computing max bottom value for
! wrho_bt_submeso. Default diag_step=1200.
! </DATA>
!
! <DATA NAME="submeso_skew_flux" TYPE="logical">
! For computing the tendency as convergence of skew flux.
! This is the recommended method.
! Default submeso_skew_flux=.true.
! </DATA>
!
! <DATA NAME="submeso_advect_flux" TYPE="logical">
! For computing the tendency as convergence of advective flux.
! This approach uses either a flux limited sweby advection or
! first order upwind, both of which ensure that the resulting
! tendency will not create extrema in the tracer field.
! Default submeso_advect_flux=.false.
! </DATA>
! <DATA NAME="submeso_advect_upwind" TYPE="logical">
! For computing the tendency as convergence of a first order
! advective flux.
! Default submeso_advect_upwind=.true.
! </DATA>
! <DATA NAME="submeso_advect_sweby" TYPE="logical">
! For computing the tendency as convergence of a sweby
! advective flux. This routine is incomplete and has a bug.
! Default submeso_advect_sweby=.false.
! </DATA>
! <DATA NAME="submeso_advect_limit" TYPE="logical">
! For limiting the value of the horizontal transports
! to be less than a velocity scale set by limit_psi_velocity_scale.
! This option is not needed if limit_psi=.true.
! Default submeso_advect_limit=.false.
! </DATA>
! <DATA NAME="submeso_advect_zero_bdy" TYPE="logical">
! For removing the advective transport next to boundaries.
! This is useful since computation of the advective transport
! velocity components can be problematic next to boundaries.
! Default submeso_advect_zero_bdy=.false.
! </DATA>
! <DATA NAME="smooth_advect_transport" TYPE="logical">
! For doing a horizontal 1-2-1 smoothing on the diagnosed
! uhrho_et_submeso and vhrho_nt_submeso fields.
! Default smooth_advect_transport=.true.
! </DATA>
! <DATA NAME="smooth_advect_transport_num" TYPE="integer">
! Number of iterations for the smooothing of horizontal transport.
! Default smooth_advect_transport_num=2.
! </DATA>
!
! <DATA NAME="submeso_diffusion" TYPE="logical">
! For computing a horizontal diffusive flux in the boundary layer
! as determined by the strength of the vector streamfunction.
! Default submeso_diffusion=.false.
! </DATA>
! <DATA NAME="submeso_diffusion_biharmonic" TYPE="logical">
! The default submeso diffusion is Laplacian. However, one may wish to
! use a biharmonic mixing operator instead.
! Default submeso_diffusion_biharmonic=.false.
! </DATA>
! <DATA NAME="submeso_diffusion_scale" UNITS="dimensionless" TYPE="real">
! A dimensionless scaling to be used for scaling up or down the effects from
! horizontal diffusion in the boundary layer. Default submeso_diffusion_scale=1.0.
! </DATA>
!
! <DATA NAME="use_hblt_constant" TYPE="logical">
! For running with a constant boundary layer depth. This for the case when
! not using a realistic mixed layer scheme. Default use_hblt_constant=.false.
! </DATA>
! <DATA NAME="constant_hblt" UNITS="metre" TYPE="real">
! The boundary layer depth for the case when use_hblt_constant=.true.
! Default constant_hblt=100.0.
! </DATA>
! <DATA NAME="use_hblt_equal_mld" TYPE="logical">
! For using the diagnosed mld as the hblt for submeso.
! This is useful for those test models that do not have a mixed layer
! scheme enabled, such as KPP, where the mixed layer scheme provides a
! boundary layer depth. In this case, it is sensible to employ the diagnosed
! mixed layer depth for the submeso scheme. Additionally, in general it is
! more physical to use the mld than the KPP hblt as the depth over which
! the submesoscale eddies act. Hence, default use_hblt_equal_mld=.true.
! </DATA>
! <DATA NAME="min_kblt" UNITS="dimensionless" TYPE="integer">
! The minimum number of vertical cells in the surface boundary layer
! that are required in order to compute the submesoscale streamfunction.
! Default min_kblt=4. Need at least three to fit a parabola with zero
! streamfunction at the top and bottom of the boundary layer.
! </DATA>
!
! <DATA NAME="minimum_hblt" TYPE="real" UNITS="metre">
! For setting a floor to the hblt used for submesoscale scheme.
! Default minimum_hblt=0.0.
! </DATA>
!
! <DATA NAME="smooth_hblt" TYPE="logical">
! For smoothing on the submeso bldepth field. This is useful
! since the bldepth obtained from KPP or diagnosed mld can
! have some grid noise.
! Default smooth_hblt=.false. since this agrees with legacy.
! Note that this scheme fails to reproduce across
! processor layout, so it remains broken.
! </DATA>
! <DATA NAME="smooth_hblt_num" TYPE="integer">
! Number of iterations for the smooothing of bldepth.
! Default smooth_hblt_num=1.
! </DATA>
!
! <DATA NAME="use_psi_legacy" TYPE="logical">
! For computing psi using older legacy methods.
! These methods are not ideal, and can be problematic
! depending on nml settings for the limiters and smoothers.
! This option is retained only for legacy purposes.
! Default use_psi_legacy=.false.
! </DATA>
! <DATA NAME="smooth_psi" TYPE="logical">
! For doing a horizontal 1-2-1 smoothing on the
! psix_horz and psiy_horz fields.
! Default smooth_psi=.true.
! </DATA>
! <DATA NAME="smooth_psi_num" TYPE="integer">
! Number of iterations for the smooothing of psi.
! Default smooth_psi_num=2.
! </DATA>
! <DATA NAME="limit_psi" TYPE="logical">
! For limiting the magnitude of psi in order to reduce possibility of
! model crashes. Rescales the full psi to maintain vertical structure
! but to keep overall magnitude within bounds.
! Default limit_psi=.false.
! </DATA>
! <DATA NAME="limit_psi_velocity_scale" UNITS="metre/sec" TYPE="real">
! Velocity scale used to limit the value of psi when limit_psi=.true.
! Default limit_psi_velocity_scale=5.0
! </DATA>
!
! <DATA NAME="submeso_limit_flux" TYPE="logical">
! For limiting the fluxes arising from submeso scheme, according to
! tmask_limit. When reach a point where tmask_limit=1.0, then set
! the submeso flux for this cell to zero.
! Default submeso_limit_flux=.true.
! </DATA>
!
! <DATA NAME="coefficient_ce" UNITS="dimensionless" TYPE="real">
! The dimensionless coefficient from the Fox-Kemper etal scheme.
! They recommend setting coefficient_ce between 0.06 and 0.08.
! Default coefficient_ce=0.07.
! </DATA>
! <DATA NAME="time_constant" UNITS="seconds" TYPE="real">
! Timescale to mix momentum across the mixed layer.
! Default time_constant=86400.0 = 1day.
! </DATA>
! <DATA NAME="front_length_const" UNITS="metre" TYPE="real">
! Take constant horizontal length scale of submesoscale front.
! Default front_length_const=5e3.
! </DATA>
! <DATA NAME="front_length_deform_radius" TYPE="logical">
! To compute the front length using the mixed layer deformation
! radius. Default front_length_deform_radius=.true. Note,
! will have a floor on the variable front length set by the
! nml setting for front_length_const.
! </DATA>
!
!</NAMELIST>
use constants_mod, only: epsln
use diag_manager_mod, only: register_diag_field, register_static_field, need_data, send_data
use fms_mod, only: write_version_number, open_namelist_file, close_file, check_nml_error
use fms_mod, only: stdout, stdlog, read_data, NOTE, FATAL, WARNING
use mpp_domains_mod, only: mpp_update_domains, XUPDATE, YUPDATE, CGRID_NE, EDGEUPDATE, SUPDATE, WUPDATE
use mpp_domains_mod, only: mpp_global_sum, NON_BITWISE_EXACT_SUM
use mpp_mod, only: input_nml_file, mpp_error, mpp_max, mpp_pe
use time_manager_mod, only: time_type, increment_time
use ocean_domains_mod, only: get_local_indices, set_ocean_domain
use ocean_operators_mod, only: FMX, FMY, FDX_T, FDY_T, BDX_ET, BDY_NT
use ocean_parameters_mod, only: missing_value, DEPTH_BASED, omega_earth, grav
use ocean_parameters_mod, only: rho0, rho0r, onehalf, onefourth, onesixth, oneeigth
use ocean_tracer_diag_mod,only: calc_mixed_layer_depth, diagnose_eta_tend_3dflux
use ocean_types_mod, only: tracer_2d_type, tracer_3d_0_nk_type, tracer_3d_1_nk_type
use ocean_types_mod, only: ocean_time_type, ocean_domain_type, ocean_grid_type, ocean_options_type
use ocean_types_mod, only: ocean_prog_tracer_type, ocean_thickness_type, ocean_density_type
use ocean_util_mod, only: write_timestamp, diagnose_2d, diagnose_3d, diagnose_sum, write_chksum_3d
use ocean_tracer_util_mod,only: diagnose_3d_rho
use ocean_workspace_mod, only: wrk1, wrk2, wrk3, wrk4, wrk5, wrk1_2d, wrk2_2d, wrk1_v
implicit none
private
! for diagnostics
real, dimension(:,:,:), allocatable :: advect_tendency
! internally set for computing watermass diagnstics
logical :: compute_watermass_diag = .false.
logical :: compute_watermass_diff_diag = .false.
! for diagnostics
integer :: id_kblt_submeso =-1
integer :: id_hblt_submeso =-1
integer :: id_mu_submeso =-1
integer :: id_dmu_submeso =-1
integer :: id_psix_submeso =-1
integer :: id_psiy_submeso =-1
integer :: id_tx_trans_submeso =-1
integer :: id_ty_trans_submeso =-1
integer :: id_tz_trans_submeso =-1
integer :: id_tx_trans_submeso_adv =-1
integer :: id_ty_trans_submeso_adv =-1
integer :: id_tz_trans_submeso_adv =-1
integer :: id_tx_trans_nrho_submeso=-1
integer :: id_ty_trans_nrho_submeso=-1
integer :: id_tz_trans_nrho_submeso=-1
integer :: id_tx_trans_rho_submeso=-1
integer :: id_ty_trans_rho_submeso=-1
integer :: id_tz_trans_rho_submeso=-1
integer :: id_front_length_submeso =-1
integer :: id_buoy_freq_ave_submeso=-1
integer :: id_uhrho_et_submeso =-1
integer :: id_vhrho_nt_submeso =-1
integer :: id_wrho_bt_submeso =-1
integer :: id_u_et_submeso =-1
integer :: id_v_nt_submeso =-1
integer :: id_w_bt_submeso =-1
integer :: id_subdiff_diffusivity =-1
integer :: id_eta_tend_submeso_flx =-1
integer :: id_eta_tend_submeso_flx_glob =-1
integer :: id_eta_tend_submeso_tend =-1
integer :: id_eta_tend_submeso_tend_glob=-1
integer :: id_eta_tend_subdiff_flx =-1
integer :: id_eta_tend_subdiff_flx_glob =-1
integer :: id_eta_tend_subdiff_tend =-1
integer :: id_eta_tend_subdiff_tend_glob=-1
integer :: id_neut_rho_submeso =-1
integer :: id_pot_rho_submeso =-1
integer :: id_neut_rho_submeso_on_nrho =-1
integer :: id_wdian_rho_submeso =-1
integer :: id_wdian_rho_submeso_on_nrho =-1
integer :: id_tform_rho_submeso =-1
integer :: id_tform_rho_submeso_on_nrho =-1
integer :: id_neut_temp_submeso =-1
integer :: id_neut_temp_submeso_on_nrho =-1
integer :: id_wdian_temp_submeso =-1
integer :: id_wdian_temp_submeso_on_nrho =-1
integer :: id_tform_temp_submeso =-1
integer :: id_tform_temp_submeso_on_nrho =-1
integer :: id_neut_salt_submeso =-1
integer :: id_neut_salt_submeso_on_nrho =-1
integer :: id_wdian_salt_submeso =-1
integer :: id_wdian_salt_submeso_on_nrho =-1
integer :: id_tform_salt_submeso =-1
integer :: id_tform_salt_submeso_on_nrho =-1
integer :: id_neut_rho_subdiff =-1
integer :: id_neut_rho_subdiff_on_nrho =-1
integer :: id_wdian_rho_subdiff =-1
integer :: id_wdian_rho_subdiff_on_nrho =-1
integer :: id_tform_rho_subdiff =-1
integer :: id_tform_rho_subdiff_on_nrho =-1
integer :: id_neut_temp_subdiff =-1
integer :: id_neut_temp_subdiff_on_nrho =-1
integer :: id_wdian_temp_subdiff =-1
integer :: id_wdian_temp_subdiff_on_nrho =-1
integer :: id_tform_temp_subdiff =-1
integer :: id_tform_temp_subdiff_on_nrho =-1
integer :: id_neut_salt_subdiff =-1
integer :: id_neut_salt_subdiff_on_nrho =-1
integer :: id_wdian_salt_subdiff =-1
integer :: id_wdian_salt_subdiff_on_nrho =-1
integer :: id_tform_salt_subdiff =-1
integer :: id_tform_salt_subdiff_on_nrho =-1
integer, dimension(:), allocatable :: id_xflux_submeso ! i-directed flux
integer, dimension(:), allocatable :: id_yflux_submeso ! j-directed flux
integer, dimension(:), allocatable :: id_zflux_submeso ! k-directed flux
integer, dimension(:), allocatable :: id_xflux_submeso_int_z ! vertically integrated i-flux
integer, dimension(:), allocatable :: id_yflux_submeso_int_z ! vertically integrated j-flux
integer, dimension(:), allocatable :: id_submeso ! tendency from streamfunction portion of submesoscale param
integer, dimension(:), allocatable :: id_submeso_on_nrho ! tendency from streamfunction portion of submesoscale param binned to neutral density
integer, dimension(:), allocatable :: id_xflux_subdiff ! i-directed horz diffusive flux
integer, dimension(:), allocatable :: id_yflux_subdiff ! j-directed horz diffusive flux
integer, dimension(:), allocatable :: id_xflux_subdiff_int_z ! vertically integrated horz diffusive i-flux
integer, dimension(:), allocatable :: id_yflux_subdiff_int_z ! vertically integrated horz diffusive j-flux
integer, dimension(:), allocatable :: id_subdiff ! tendency from horz diffusion associated with submesoscale
logical :: used
#include <ocean_memory.h>
type(ocean_domain_type), pointer :: Dom => NULL()
type(ocean_grid_type), pointer :: Grd => NULL()
type(ocean_domain_type), save :: Dom_flux_sub
character(len=128) :: version='$$'
character (len=128) :: tagname = '$Name: tikal $'
#ifdef MOM_STATIC_ARRAYS
real, dimension(isd:ied,jsd:jed,nk) :: uhrho_et_submeso ! i-component of transport for submeso
real, dimension(isd:ied,jsd:jed,nk) :: vhrho_nt_submeso ! j-component of transport for submeso
real, dimension(isd:ied,jsd:jed,0:nk) :: wrho_bt_submeso ! k-component of transport for submeso
real, dimension(isd:ied,jsd:jed,nk,0:1) :: delqc !density weighted (kg/m^3) quarter cell thickness(m)
real, dimension(isd:ied,jsd:jed,0:nk) :: dzwtr !(1/dzwt)(m^-1)
real, dimension(isd:ied,jsd:jed,0:1) :: dtew !grid distances from T-point to cell faces (m)
real, dimension(isd:ied,jsd:jed,0:1) :: dtns !grid distances from T-point to cell faces (m)
real, dimension(isd:ied,jsd:jed,0:1) :: dtwedyt !horizontal areas (m^2) of quarter cell
real, dimension(isd:ied,jsd:jed,0:1) :: dxtdtsn !horizontal areas (m^2) of quarter cell
real, dimension(isd:ied,jsd:jed,nk,0:1) :: psix ! streamfunction x-component (m^2/sec)
real, dimension(isd:ied,jsd:jed,nk,0:1) :: psiy ! streamfunction y-component (m^2/sec)
real, dimension(isd:ied,jsd:jed) :: hblt ! boundary layer depth (m)
real, dimension(isd:ied,jsd:jed) :: grid_length ! grid length scale (m)
real, dimension(isd:ied,jsd:jed) :: front_length_inv ! inverse front length (1/m)
real, dimension(isd:ied,jsd:jed) :: buoy_freq_ave ! buoyancy frequency averaged over mixed layer depth (1/sec)
real, dimension(isd:ied,jsd:jed) :: time_factor ! time factor (sec) for computing streamfunction
real, dimension(isd:ied,jsd:jed) :: coriolis_param ! absolute value of the Coriolis parameter (sec^-1) on T-grid
integer, dimension(isd:ied,jsd:jed) :: kblt ! k-level encompassing hblt
real, dimension(isd:ied,jsd:jed,nk) :: flux_x ! i-component to tracer flux
real, dimension(isd:ied,jsd:jed,nk) :: flux_y ! j-component to tracer flux
real, dimension(isd:ied,jsd:jed,0:nk) :: flux_z ! k-component to tracer flux
#else
real, dimension(:,:,:), allocatable :: uhrho_et_submeso ! i-component of transport for submeso
real, dimension(:,:,:), allocatable :: vhrho_nt_submeso ! j-component of transport for submeso
real, dimension(:,:,:), allocatable :: wrho_bt_submeso ! k-component of transport for submeso
real, dimension(:,:,:,:), allocatable :: delqc !density weighted (kg/m^3) quarter cell thickness(m)
real, dimension(:,:,:), allocatable :: dzwtr !(1/dzwt)(m^-1)
real, dimension(:,:,:), allocatable :: dtew !grid distances from T-point to cell faces (m)
real, dimension(:,:,:), allocatable :: dtns !grid distances from T-point to cell faces (m)
real, dimension(:,:,:), allocatable :: dtwedyt !horizontal areas (m^2) of quarter cell
real, dimension(:,:,:), allocatable :: dxtdtsn !horizontal areas (m^2) of quarter cell
real, dimension(:,:,:,:), allocatable :: psix ! streamfunction x-component (m^2/sec)
real, dimension(:,:,:,:), allocatable :: psiy ! streamfunction y-component (m^2/sec)
real, dimension(:,:), allocatable :: hblt ! boundary layer depth (m)
real, dimension(:,:), allocatable :: grid_length ! grid length scale (m)
real, dimension(:,:), allocatable :: front_length_inv ! inverse front length (1/m)
real, dimension(:,:), allocatable :: buoy_freq_ave ! buoyancy frequency averaged over mixed layer depth (1/sec)
real, dimension(:,:), allocatable :: time_factor ! time factor (sec) for computing streamfunction
real, dimension(:,:), allocatable :: coriolis_param ! absolute value of the Coriolis parameter (sec^-1) on T-grid
integer, dimension(:,:), allocatable :: kblt ! k-level encompassing hblt
real, dimension(:,:,:), allocatable :: flux_x ! i-component to tracer flux
real, dimension(:,:,:), allocatable :: flux_y ! j-component to tracer flux
real, dimension(:,:,:), allocatable :: flux_z ! k-component to tracer flux
#endif
! for saving the horizontally dependent portion of psix and psiy
real, dimension(:,:,:), allocatable :: psix_horz ! streamfunction x-component sans vertical structure (m^2/sec)
real, dimension(:,:,:), allocatable :: psiy_horz ! streamfunction y-component sans vertical structure (m^2/sec)
! for eta_tend diagnostics
real, dimension(:,:,:), allocatable :: flux_x_temp
real, dimension(:,:,:), allocatable :: flux_y_temp
real, dimension(:,:,:), allocatable :: flux_z_temp
real, dimension(:,:,:), allocatable :: flux_x_salt
real, dimension(:,:,:), allocatable :: flux_y_salt
real, dimension(:,:,:), allocatable :: flux_z_salt
! for advecting tracers with sweby submesoscale advection
real, dimension(:,:,:), allocatable :: tmask_mdfl
real, dimension(:,:,:), allocatable :: tracer_mdfl
type :: tracer_mdfl_type
real, dimension(:,:,:), pointer :: field => NULL()
end type tracer_mdfl_type
type(tracer_mdfl_type), dimension(:), allocatable :: tracer_mdfl_all ! tracer array for sweby advection
type(ocean_domain_type), save :: Dom_mdfl
type(tracer_3d_1_nk_type), dimension(:), allocatable :: dTdx ! tracer partial derivative (tracer/m)
type(tracer_3d_1_nk_type), dimension(:), allocatable :: dTdy ! tracer partial derivative (tracer/m)
type(tracer_3d_0_nk_type), dimension(:), allocatable :: dTdz ! tracer partial derivative (tracer/m)
type(tracer_3d_1_nk_type), save :: dSdx ! density-salinity partial derivative (tracer/m)
type(tracer_3d_1_nk_type), save :: dSdy ! density-salinity partial derivative (tracer/m)
type(tracer_3d_0_nk_type), save :: dSdz ! density-salinity partial derivative (tracer/m)
public ocean_submesoscale_init
public submeso_restrat
private tracer_derivs
private salinity_derivs
private compute_flux_x
private compute_flux_y
private compute_flux_z
private compute_psi
private compute_psi_legacy
private compute_bldepth
private compute_advect_transport
private compute_submeso_skewsion
private compute_submeso_upwind
private compute_submeso_sweby
private maximum_bottom_w_submeso
private transport_on_rho_submeso
private transport_on_nrho_submeso
private transport_on_nrho_submeso_adv
private watermass_diag_init
private watermass_diag
private watermass_diag_diffusion
integer :: index_temp=-1
integer :: index_salt=-1
integer :: num_prog_tracers=0
! for diagnosing fluxes
real :: flux_sign
! vertical coordinate
integer :: vert_coordinate_class=1
! for output
integer :: unit=6
real :: fiveover21
real :: eightover21
real :: dtime
real :: ce_grav_rho0r
real :: time_constant2_r
real :: grav_rho0r
real :: front_length_const_inv
real :: front_length_max=1e10
! for global area normalization
real :: cellarea_r
logical :: module_is_initialized=.false.
! for specifying transport units
! can either be Sv or mks
character(len=32) :: transport_dims ='Sv (10^9 kg/s)'
real :: transport_convert=1.0e-9
! for eta_tend diagnostics (internally set)
logical :: diagnose_eta_tend_submeso_flx=.false.
logical :: diagnose_eta_tend_subdiff_flx=.false.
! internally determined
logical :: diag_advect_transport = .false.
! nml parameters
logical :: use_this_module = .false.
logical :: debug_this_module = .false.
logical :: submeso_skew_flux = .true.
logical :: submeso_advect_flux = .false.
logical :: submeso_advect_upwind = .true.
logical :: submeso_advect_sweby = .false.
logical :: submeso_advect_limit = .false.
logical :: submeso_advect_zero_bdy = .false.
logical :: submeso_limit_flux = .true.
logical :: use_hblt_constant = .false.
logical :: use_hblt_equal_mld = .true.
logical :: smooth_hblt = .false.
logical :: smooth_psi = .true.
logical :: smooth_advect_transport = .true.
logical :: front_length_deform_radius = .true.
logical :: limit_psi = .true.
logical :: use_psi_legacy = .false.
logical :: submeso_diffusion = .false.
logical :: submeso_diffusion_biharmonic = .false.
real :: limit_psi_velocity_scale = 0.5
real :: coefficient_ce = 0.07
real :: time_constant = 86400.0
real :: front_length_const = 5e3
real :: constant_hblt = 100.0
real :: minimum_hblt = 0.0
real :: submeso_diffusion_scale = 1.0
integer :: min_kblt = 4
integer :: diag_step = 1200
integer :: smooth_psi_num = 2
integer :: smooth_advect_transport_num = 2
integer :: smooth_hblt_num = 2
namelist /ocean_submesoscale_nml/ use_this_module, debug_this_module, diag_step, &
use_hblt_constant, use_hblt_equal_mld, &
smooth_hblt, smooth_hblt_num, constant_hblt, &
coefficient_ce, time_constant, front_length_const, &
min_kblt, minimum_hblt, smooth_psi, smooth_psi_num, &
front_length_deform_radius, &
limit_psi, use_psi_legacy, limit_psi_velocity_scale, &
submeso_limit_flux, smooth_advect_transport, &
smooth_advect_transport_num, &
submeso_skew_flux, submeso_advect_flux, &
submeso_advect_upwind, submeso_advect_sweby, &
submeso_advect_limit, submeso_advect_zero_bdy, &
submeso_diffusion, submeso_diffusion_scale, &
submeso_diffusion_biharmonic
contains
!#######################################################################
! <SUBROUTINE NAME="ocean_submesoscale_init">
!
! <DESCRIPTION>
! Initialization for the ocean_submesoscale module.
! </DESCRIPTION>
subroutine ocean_submesoscale_init(Grid, Domain, Time, Dens, T_prog, Ocean_options, dtime_t, &
ver_coordinate_class, cmip_units, debug)
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_density_type), intent(in) :: Dens
type(ocean_prog_tracer_type), intent(in) :: T_prog(:)
type(ocean_options_type), intent(inout) :: Ocean_options
real, intent(in) :: dtime_t
integer, intent(in) :: ver_coordinate_class
logical, intent(in) :: cmip_units
logical, intent(in), optional :: debug
integer :: unit, io_status, ierr
integer :: i,j,k,n
integer :: num_methods=0
integer :: stdoutunit,stdlogunit
stdoutunit=stdout();stdlogunit=stdlog()
if ( module_is_initialized ) then
call mpp_error(FATAL,&
'==>Error from ocean_submesoscale_init_mod: module already initialized')
endif
module_is_initialized = .TRUE.
call write_version_number(version, tagname)
#ifdef INTERNAL_FILE_NML
read (input_nml_file, nml=ocean_submesoscale_nml, iostat=io_status)
ierr = check_nml_error(io_status,'ocean_submesoscale_nml')
#else
unit = open_namelist_file()
read(unit, ocean_submesoscale_nml,iostat=io_status)
ierr = check_nml_error(io_status, 'ocean_submesoscale_nml')
call close_file(unit)
#endif
write (stdoutunit,'(/)')
write (stdoutunit,ocean_submesoscale_nml)
write (stdlogunit,ocean_submesoscale_nml)
Dom => Domain
Grd => Grid
#ifndef MOM_STATIC_ARRAYS
call get_local_indices(Domain, isd, ied, jsd, jed, isc, iec, jsc, jec)
nk = Grid%nk
#endif
if(cmip_units) then
transport_convert=1.0
transport_dims = 'kg/s'
else
transport_convert=1.0e-9
transport_dims = 'Sv (10^9 kg/s)'
endif
if (PRESENT(debug) .and. .not. debug_this_module) then
debug_this_module = debug
endif
if(debug_this_module) then
write(stdoutunit,'(a)') '==>Note: running ocean_submesoscale_mod with debug_this_module=.true.'
endif
if(use_this_module) then
call mpp_error(NOTE, '==>Note: USING ocean_submesoscale_mod')
Ocean_options%submesoscale = 'Used submesoscale closure for surface restratification.'
else
call mpp_error(NOTE, '==>Note: NOT using ocean_submesoscale_mod')
Ocean_options%submesoscale = 'Did NOT use submesoscale closure for surface restratification.'
return
endif
if(use_hblt_equal_mld) then
write(stdoutunit,'(a)') &
'==>Note: For ocean_submesoscale, setting bldepth equal to diagnosed mld.'
endif
if(use_hblt_constant) then
write(stdoutunit,'(a)') &
'==>Note: For ocean_submesoscale, setting bldepth equal to prescribed constant.'
endif
if(use_hblt_equal_mld .and. use_hblt_constant) then
call mpp_error(FATAL, &
'==>Error: in ocean_submesoscale, use_hblt_equal_mld & use_hblt_constant cannot both be true.')
endif
if(submeso_advect_flux) then
num_methods=num_methods+1
flux_sign = 1.0
if(submeso_advect_upwind) then
write(stdoutunit,'(a)') &
'==>Note: For ocean_submesoscale, tendency computed as convergence of upwind advective flux.'
elseif(submeso_advect_sweby) then
write(stdoutunit,'(a)') &
'==>Note: For ocean_submesoscale, tendency computed as convergence of mdfl_sweby advective flux; has known bugs!'
else
call mpp_error(FATAL, &
'==>Error: in ocean_submesoscale: submeso_advect_flux=.true. yet no advection scheme chosen')
endif
if(use_psi_legacy) then
call mpp_error(FATAL, &
'==>Error: in ocean_submesoscale: submeso_advect_flux=.true. is not compatible with use_psi_legacy=.true.')
endif
endif
if(submeso_skew_flux) then
write(stdoutunit,'(a)') &
'==>Note: For ocean_submesoscale, tendency computed as skew flux convergence.'
num_methods=num_methods+1
flux_sign = -1.0
endif
if(num_methods > 1) then
call mpp_error(FATAL, &
'==>Error: in ocean_submesoscale, can choose only one method for computing tendency.')
endif
if(num_methods == 0) then
call mpp_error(FATAL, &
'==>Error: in ocean_submesoscale, must choose one method for computing tendency.')
endif
if(use_psi_legacy) then
write(stdoutunit,'(a)') &
'==>Note: in ocean_submesoscale: use_psi_legacy=.true. This method has problems, and is retained solely for legacy.'
endif
if(submeso_diffusion) then
if(submeso_diffusion_biharmonic) then
write(stdoutunit,'(a)') &
'==>Note: in ocean_submesoscale: submeso_diffusion=.true. Adding horizontal biharmonic mixing in blayer.'
else
write(stdoutunit,'(a)') &
'==>Note: in ocean_submesoscale: submeso_diffusion=.true. Adding horizontal Laplacian diffusion in blayer.'
endif
endif
fiveover21 = 5.0/21.0
eightover21 = 8.0/21.0
grav_rho0r = grav*rho0r
ce_grav_rho0r = coefficient_ce*grav*rho0r
time_constant2_r = 1.0/(time_constant**2)
dtime = dtime_t
vert_coordinate_class = ver_coordinate_class
front_length_const_inv = 1.0/front_length_const
cellarea_r = 1.0/(epsln + Grd%tcellsurf)
call set_ocean_domain(Dom_mdfl,Grd,xhalo=2,yhalo=2,name='mdfl',maskmap=Dom%maskmap)
call set_ocean_domain(Dom_flux_sub, Grd,xhalo=Dom%xhalo,yhalo=Dom%yhalo,name='flux dom submeso',&
maskmap=Dom%maskmap)
! for diagnostics
allocate( advect_tendency(isd:ied,jsd:jed,nk) )
advect_tendency(:,:,:) = 0.0
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( dTdx(num_prog_tracers) )
allocate( dTdy(num_prog_tracers) )
allocate( dTdz(num_prog_tracers) )
#ifndef MOM_STATIC_ARRAYS
allocate (uhrho_et_submeso(isd:ied,jsd:jed,nk))
allocate (vhrho_nt_submeso(isd:ied,jsd:jed,nk))
allocate (wrho_bt_submeso(isd:ied,jsd:jed,0:nk))
allocate (delqc(isd:ied,jsd:jed,nk,0:1))
allocate (dzwtr(isd:ied,jsd:jed,0:nk))
allocate (dtew(isd:ied,jsd:jed,0:1))
allocate (dtns(isd:ied,jsd:jed,0:1))
allocate (dtwedyt(isd:ied,jsd:jed,0:1))
allocate (dxtdtsn(isd:ied,jsd:jed,0:1))
allocate (psix(isd:ied,jsd:jed,nk,0:1))
allocate (psiy(isd:ied,jsd:jed,nk,0:1))
allocate (hblt(isd:ied,jsd:jed))
allocate (grid_length(isd:ied,jsd:jed))
allocate (front_length_inv(isd:ied,jsd:jed))
allocate (buoy_freq_ave(isd:ied,jsd:jed))
allocate (time_factor(isd:ied,jsd:jed))
allocate (coriolis_param(isd:ied,jsd:jed))
allocate (kblt(isd:ied,jsd:jed))
allocate (flux_x(isd:ied,jsd:jed,nk) )
allocate (flux_y(isd:ied,jsd:jed,nk) )
allocate (flux_z(isd:ied,jsd:jed,0:nk) )
do n=1,num_prog_tracers
allocate ( dTdx(n)%field(isd:ied,jsd:jed,nk) )
allocate ( dTdy(n)%field(isd:ied,jsd:jed,nk) )
allocate ( dTdz(n)%field(isd:ied,jsd:jed,0:nk) )
enddo
allocate ( dSdx%field(isd:ied,jsd:jed,nk) )
allocate ( dSdy%field(isd:ied,jsd:jed,nk) )
allocate ( dSdz%field(isd:ied,jsd:jed,0:nk) )
#endif
allocate (psix_horz(isd:ied,jsd:jed,0:1))
allocate (psiy_horz(isd:ied,jsd:jed,0:1))
uhrho_et_submeso = 0.0
vhrho_nt_submeso = 0.0
wrho_bt_submeso = 0.0
do n=1,num_prog_tracers
dTdx(n)%field(:,:,:) = 0.0
dTdy(n)%field(:,:,:) = 0.0
dTdz(n)%field(:,:,:) = 0.0
enddo
dSdx%field(:,:,:) = 0.0
dSdy%field(:,:,:) = 0.0
dSdz%field(:,:,:) = 0.0
psix = 0.0
psiy = 0.0
psix_horz = 0.0
psiy_horz = 0.0
kblt = 0
hblt = 0.0
flux_x = 0.0
flux_y = 0.0
flux_z = 0.0
dzwtr = 0.0
delqc = 0.0
grid_length(:,:) = 2.0*Grd%dxt(:,:)*Grd%dyt(:,:)/(Grd%dxt(:,:)+Grd%dyt(:,:))
front_length_inv(:,:) = front_length_const_inv
buoy_freq_ave(:,:) = 0.0
coriolis_param(:,:) = 2.0*omega_earth*abs(sin(Grd%phit(:,:)))
dtew(:,:,0) = Grd%dtw(:,:)
dtew(:,:,1) = Grd%dte(:,:)
dtns(:,:,0) = Grd%dts(:,:)
dtns(:,:,1) = Grd%dtn(:,:)
dtwedyt(:,:,:) = 0.0
dtwedyt(:,:,0) = Grd%dte(:,:)*Grd%dyt(:,:)
do i=isc-1,iec
dtwedyt(i,:,1) = Grd%dtw(i+1,:)*Grd%dyt(i+1,:)
enddo
dxtdtsn(:,:,:) = 0.0
dxtdtsn(:,:,0) = Grd%dxt(:,:)*Grd%dtn(:,:)
do j=jsc-1,jec
dxtdtsn(:,j,1) = Grd%dxt(:,j+1)*Grd%dts(:,j+1)
enddo
do j=jsd,jed
do i=isd,ied
time_factor(i,j) = 1.0/(sqrt( time_constant2_r + coriolis_param(i,j)**2 ))
enddo
enddo
! for the sweby advection scheme
if(submeso_advect_flux .and. submeso_advect_sweby) then
allocate(tracer_mdfl_all(num_prog_tracers))
allocate(tmask_mdfl (isc-2:iec+2,jsc-2:jec+2,nk))
allocate(tracer_mdfl(isc-2:iec+2,jsc-2:jec+2,nk))
do n=1,num_prog_tracers
allocate (tracer_mdfl_all(n)%field(isc-2:iec+2,jsc-2:jec+2,nk))
enddo
tmask_mdfl = 0.0
tracer_mdfl = 0.0
do n=1,num_prog_tracers
tracer_mdfl_all(n)%field(:,:,:) = 0.0
enddo
do k=1,nk
do j=jsc,jec
do i=isc,iec
tmask_mdfl(i,j,k) = Grd%tmask(i,j,k)
enddo
enddo
enddo
call mpp_update_domains(tmask_mdfl,Dom_mdfl%domain2d)
endif
! diagnostics
id_kblt_submeso = register_diag_field ('ocean_model', 'kblt_submeso', &
Grid%tracer_axes(1:2), Time%model_time, &
'Number of k-levels in boundary layer for submesoscale closure', &
'dimensionless', missing_value=missing_value, range=(/-1.0,1e6/))
id_hblt_submeso = register_diag_field ('ocean_model', 'hblt_submeso', &
Grid%tracer_axes(1:2), Time%model_time, &
'Boundary layer depth used for submesoscale closure', &
'metre', missing_value=missing_value, range=(/-1.0,1e6/))
id_front_length_submeso = register_diag_field ('ocean_model', 'front_length_submeso', &
Grid%tracer_axes(1:2), Time%model_time, &
'Front length used for submesoscale closure', &
'metre', missing_value=missing_value, range=(/-1.0,1e6/))
id_buoy_freq_ave_submeso = register_diag_field ('ocean_model', 'buoy_freq_ave_submeso',&
Grid%tracer_axes(1:2), Time%model_time, &
'Buoyancy frequency averaged over depth of mixed layer for submesoscale closure', &
'1/sec', missing_value=missing_value, range=(/-1.0,1e6/))
id_mu_submeso = register_diag_field ('ocean_model', 'mu_submeso', &
Grd%tracer_axes(1:3), Time%model_time, &
'vertical structure function for submesoscale streamfunction', &
'dimensionless', missing_value=missing_value, range=(/-1.e2,1e2/))
id_dmu_submeso = register_diag_field ('ocean_model', 'dmu_submeso', &
Grd%tracer_axes(1:3), Time%model_time, &
'vertical derivative of vertical structure function for submesoscale streamfunction', &
'1/metre', missing_value=missing_value, range=(/-1.e6,1e6/))
id_psix_submeso = register_diag_field ('ocean_model', 'psix_submeso', &
Grd%tracer_axes_flux_y(1:3), Time%model_time, &
'i-comp of submesoscale streamfunction', &
'm^2/sec', missing_value=missing_value, range=(/-1.e10,1e10/))
id_psiy_submeso = register_diag_field ('ocean_model', 'psiy_submeso', &
Grd%tracer_axes_flux_x(1:3), Time%model_time, &
'j-comp of submesoscale streamfunction', &
'm^2/sec', missing_value=missing_value, range=(/-1.e10,1e10/))
id_tx_trans_submeso = register_diag_field ('ocean_model', 'tx_trans_submeso', &
Grd%tracer_axes_flux_x(1:3), Time%model_time, &
'T-cell mass i-transport from submesoscale param', &
trim(transport_dims), missing_value=missing_value, range=(/-1.e20,1e20/))
id_ty_trans_submeso = register_diag_field ('ocean_model', 'ty_trans_submeso', &
Grd%tracer_axes_flux_y(1:3), Time%model_time, &
'T-cell mass j-transport from submesoscale param', &
trim(transport_dims), missing_value=missing_value, range=(/-1.e20,1e20/))
id_tz_trans_submeso = register_diag_field ('ocean_model', 'tz_trans_submeso', &
Grd%tracer_axes_wt(1:3), Time%model_time, &
'T-cell mass k-transport from submesoscale param', &
trim(transport_dims), missing_value=missing_value, range=(/-1.e20,1e20/))
id_tx_trans_nrho_submeso = register_diag_field ('ocean_model','tx_trans_nrho_submeso',&
Dens%neutralrho_axes_flux_x(1:3),Time%model_time, &
'T-cell i-mass transport from submesoscale param on neutral rho' &
,trim(transport_dims),missing_value=missing_value, range=(/-1e20,1e20/))
id_ty_trans_nrho_submeso = register_diag_field ('ocean_model','ty_trans_nrho_submeso',&
Dens%neutralrho_axes_flux_y(1:3),Time%model_time, &
'T-cell j-mass transport from submesoscale param on neutral rho' &
,trim(transport_dims),missing_value=missing_value, range=(/-1e20,1e20/))
id_tz_trans_nrho_submeso = register_diag_field ('ocean_model','tz_trans_nrho_submeso',&
Dens%neutralrho_axes(1:3),Time%model_time, &
'T-cell k-mass transport from submesoscale param on neutral rho' &
,trim(transport_dims),missing_value=missing_value, range=(/-1e20,1e20/))
id_tx_trans_rho_submeso = register_diag_field('ocean_model','tx_trans_rho_submeso', &
Dens%potrho_axes_flux_x(1:3),Time%model_time, &
'T-cell i-mass transport from submesoscale param on pot_rho' &
,trim(transport_dims),missing_value=missing_value,range=(/-1e20,1e20/))
id_ty_trans_rho_submeso = register_diag_field('ocean_model','ty_trans_rho_submeso', &
Dens%potrho_axes_flux_y(1:3),Time%model_time, &
'T-cell j-mass transport from submesoscale param on pot_rho' &
,trim(transport_dims),missing_value=missing_value,range=(/-1e20,1e20/))
id_u_et_submeso = register_diag_field ('ocean_model', 'u_et_submeso', &
Grd%tracer_axes_flux_x(1:3), Time%model_time, &
'i-component of submesoscale transport velocity', &
'm/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
id_v_nt_submeso = register_diag_field ('ocean_model', 'v_nt_submeso', &
Grd%tracer_axes_flux_y(1:3), Time%model_time, &
'j-component of submesoscale transport velocity', &
'm/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
id_w_bt_submeso = register_diag_field ('ocean_model', 'w_bt_submeso', &
Grd%tracer_axes_wt(1:3), Time%model_time, &
'vertical component of submesoscale transport velocity', &
'm/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
id_uhrho_et_submeso = register_diag_field ('ocean_model', 'uhrho_et_submeso', &
Grd%tracer_axes_flux_x(1:3), Time%model_time, &
'i-component of submesoscale advective mass transport', &
'(kg/m^3)*m^2/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_uhrho_et_submeso > 0) diag_advect_transport = .true.
id_vhrho_nt_submeso = register_diag_field ('ocean_model', 'vhrho_nt_submeso', &
Grd%tracer_axes_flux_y(1:3), Time%model_time, &
'j-component of submesoscale advective mass transport', &
'(kg/m^3)*m^2/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_vhrho_nt_submeso > 0) diag_advect_transport = .true.
id_wrho_bt_submeso = register_diag_field ('ocean_model', 'wrho_bt_submeso', &
Grd%tracer_axes_wt(1:3), Time%model_time, &
'k-component of submesoscale advective mass transport', &
'(kg/m^3)*m/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_wrho_bt_submeso > 0) diag_advect_transport = .true.
id_tx_trans_submeso_adv = register_diag_field ('ocean_model', 'tx_trans_submeso_adv', &
Grd%tracer_axes_flux_x(1:3), Time%model_time, &
'i-component of submesoscale advective mass transport', &
'kg/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_tx_trans_submeso_adv > 0) diag_advect_transport = .true.
id_ty_trans_submeso_adv = register_diag_field ('ocean_model', 'ty_trans_submeso_adv', &
Grd%tracer_axes_flux_y(1:3), Time%model_time, &
'j-component of submesoscale advective mass transport', &
'kg/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_ty_trans_submeso_adv > 0) diag_advect_transport = .true.
id_tz_trans_submeso_adv = register_diag_field ('ocean_model', 'tz_trans_submeso_adv', &
Grd%tracer_axes_wt(1:3), Time%model_time, &
'k-component of submesoscale advective mass transport', &
'kg/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if(id_tz_trans_submeso_adv > 0) diag_advect_transport = .true.
if(submeso_diffusion_biharmonic) then
id_subdiff_diffusivity = register_diag_field ('ocean_model', 'subdiff_diffusivity', &
Grd%tracer_axes(1:3), Time%model_time, &
'horizontal diffusivity used for submesoscale biharmonic mixing scheme', &
'm^4/sec', missing_value=missing_value, range=(/-1.e10,1e10/))
else
id_subdiff_diffusivity = register_diag_field ('ocean_model', 'subdiff_diffusivity', &
Grd%tracer_axes(1:3), Time%model_time, &
'horizontal diffusivity used for submesoscale laplacian diffusion scheme', &
'm^2/sec', missing_value=missing_value, range=(/-1.e10,1e10/))
endif
call watermass_diag_init(Time, Dens)
id_eta_tend_submeso_flx= -1
id_eta_tend_submeso_flx= register_diag_field ('ocean_model','eta_tend_submeso_flx',&
Grd%tracer_axes(1:2), Time%model_time, &
'non-Bouss steric sea level tendency from submesoscale fluxes', 'm/s', &
missing_value=missing_value, range=(/-1e10,1.e10/))
if(id_eta_tend_submeso_flx > 0) diagnose_eta_tend_submeso_flx=.true.
id_eta_tend_submeso_flx_glob= -1
id_eta_tend_submeso_flx_glob= register_diag_field ('ocean_model', 'eta_tend_submeso_flx_glob',&
Time%model_time, &
'global mean non-bouss steric sea level tendency from submesoscale fluxes', &
'm/s', missing_value=missing_value, range=(/-1e10,1.e10/))
if(id_eta_tend_submeso_flx_glob > 0) diagnose_eta_tend_submeso_flx=.true.