MOM6
Data Types | Functions/Subroutines
mom_regularize_layers Module Reference

Detailed Description

Provides regularization of layers in isopycnal mode.

Data Types

type  regularize_layers_cs
 This control structure holds parameters used by the MOM_regularize_layers module. More...
 

Functions/Subroutines

subroutine, public regularize_layers (h, tv, dt, ea, eb, G, GV, US, CS)
 This subroutine partially steps the bulk mixed layer model. The following processes are executed, in the order listed. More...
 
subroutine regularize_surface (h, tv, dt, ea, eb, G, GV, US, CS)
 This subroutine ensures that there is a degree of horizontal smoothness in the depths of the near-surface interfaces. More...
 
subroutine find_deficit_ratios (e, def_rat_u, def_rat_v, G, GV, CS, def_rat_u_2lay, def_rat_v_2lay, halo, h)
 This subroutine determines the amount by which the harmonic mean thickness at velocity points differ from the arithmetic means, relative to the the arithmetic means, after eliminating thickness variations that are solely due to topography and aggregating all interior layers into one. More...
 
subroutine, public regularize_layers_init (Time, G, GV, param_file, diag, CS)
 Initializes the regularize_layers control structure. More...
 

Variables

integer id_clock_pass
 Clock IDs. More...
 
integer id_clock_eos
 Clock IDs. More...
 

Function/Subroutine Documentation

◆ find_deficit_ratios()

subroutine mom_regularize_layers::find_deficit_ratios ( real, dimension(szi_(g),szj_(g),szk_(g)+1), intent(in)  e,
real, dimension(szib_(g),szj_(g)), intent(out)  def_rat_u,
real, dimension(szi_(g),szjb_(g)), intent(out)  def_rat_v,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(regularize_layers_cs), pointer  CS,
real, dimension(szib_(g),szj_(g)), intent(out), optional  def_rat_u_2lay,
real, dimension(szi_(g),szjb_(g)), intent(out), optional  def_rat_v_2lay,
integer, intent(in), optional  halo,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in), optional  h 
)
private

This subroutine determines the amount by which the harmonic mean thickness at velocity points differ from the arithmetic means, relative to the the arithmetic means, after eliminating thickness variations that are solely due to topography and aggregating all interior layers into one.

Parameters
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]eInterface depths [H ~> m or kg m-2]
[out]def_rat_uThe thickness deficit ratio at u points,
[out]def_rat_vThe thickness deficit ratio at v points,
csThe control structure returned by a previous call to regularize_layers_init.
[out]def_rat_u_2layThe thickness deficit ratio at u
[out]def_rat_v_2layThe thickness deficit ratio at v
[in]haloAn extra-wide halo size, 0 by default.
[in]hLayer thicknesses [H ~> m or kg m-2].

Definition at line 617 of file MOM_regularize_layers.F90.

617  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
618  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
619  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), &
620  intent(in) :: e !< Interface depths [H ~> m or kg m-2]
621  real, dimension(SZIB_(G),SZJ_(G)), &
622  intent(out) :: def_rat_u !< The thickness deficit ratio at u points,
623  !! [nondim].
624  real, dimension(SZI_(G),SZJB_(G)), &
625  intent(out) :: def_rat_v !< The thickness deficit ratio at v points,
626  !! [nondim].
627  type(regularize_layers_cs), pointer :: cs !< The control structure returned by a
628  !! previous call to regularize_layers_init.
629  real, dimension(SZIB_(G),SZJ_(G)), &
630  optional, intent(out) :: def_rat_u_2lay !< The thickness deficit ratio at u
631  !! points when the mixed and buffer layers
632  !! are aggregated into 1 layer [nondim].
633  real, dimension(SZI_(G),SZJB_(G)), &
634  optional, intent(out) :: def_rat_v_2lay !< The thickness deficit ratio at v
635  !! pointswhen the mixed and buffer layers
636  !! are aggregated into 1 layer [nondim].
637  integer, optional, intent(in) :: halo !< An extra-wide halo size, 0 by default.
638  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
639  optional, intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
640  !! If h is not present, vertical differences
641  !! in interface heights are used instead.
642  ! Local variables
643  real, dimension(SZIB_(G),SZJ_(G)) :: &
644  h_def_u, & ! The vertically summed thickness deficits at u-points [H ~> m or kg m-2].
645  h_norm_u, & ! The vertically summed arithmetic mean thickness by which
646  ! h_def_u is normalized [H ~> m or kg m-2].
647  h_def2_u
648  real, dimension(SZI_(G),SZJB_(G)) :: &
649  h_def_v, & ! The vertically summed thickness deficits at v-points [H ~> m or kg m-2].
650  h_norm_v, & ! The vertically summed arithmetic mean thickness by which
651  ! h_def_v is normalized [H ~> m or kg m-2].
652  h_def2_v
653  real :: h_neglect ! A thickness that is so small it is usually lost
654  ! in roundoff and can be neglected [H ~> m or kg m-2].
655  real :: hmix_min ! A local copy of CS%Hmix_min [H ~> m or kg m-2].
656  real :: h1, h2 ! Temporary thicknesses [H ~> m or kg m-2].
657  integer :: i, j, k, is, ie, js, je, nz, nkmb
658 
659  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
660  if (present(halo)) then
661  is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
662  endif
663  nkmb = gv%nk_rho_varies
664  h_neglect = gv%H_subroundoff
665  hmix_min = cs%Hmix_min
666 
667  ! Determine which zonal faces are problematic.
668  do j=js,je ; do i=is-1,ie
669  ! Aggregate all water below the mixed and buffer layers for the purposes of
670  ! this diagnostic.
671  h1 = e(i,j,nkmb+1)-e(i,j,nz+1) ; h2 = e(i+1,j,nkmb+1)-e(i+1,j,nz+1)
672  if (e(i,j,nz+1) < e(i+1,j,nz+1)) then
673  if (h1 > h2) h1 = max(e(i,j,nkmb+1)-e(i+1,j,nz+1), h2)
674  elseif (e(i+1,j,nz+1) < e(i,j,nz+1)) then
675  if (h2 > h1) h2 = max(e(i+1,j,nkmb+1)-e(i,j,nz+1), h1)
676  endif
677  h_def_u(i,j) = 0.5*(h1-h2)**2 / ((h1 + h2) + h_neglect)
678  h_norm_u(i,j) = 0.5*(h1+h2)
679  enddo ; enddo
680  if (present(def_rat_u_2lay)) then ; do j=js,je ; do i=is-1,ie
681  ! This is a particular metric of the aggregation into two layers.
682  h1 = e(i,j,1)-e(i,j,nkmb+1) ; h2 = e(i+1,j,1)-e(i+1,j,nkmb+1)
683  if (e(i,j,nkmb+1) < e(i+1,j,nz+1)) then
684  if (h1 > h2) h1 = max(e(i,j,1)-e(i+1,j,nz+1), h2)
685  elseif (e(i+1,j,nkmb+1) < e(i,j,nz+1)) then
686  if (h2 > h1) h2 = max(e(i+1,j,1)-e(i,j,nz+1), h1)
687  endif
688  h_def2_u(i,j) = h_def_u(i,j) + 0.5*(h1-h2)**2 / ((h1 + h2) + h_neglect)
689  enddo ; enddo ; endif
690  do k=1,nkmb ; do j=js,je ; do i=is-1,ie
691  if (present(h)) then
692  h1 = h(i,j,k) ; h2 = h(i+1,j,k)
693  else
694  h1 = e(i,j,k)-e(i,j,k+1) ; h2 = e(i+1,j,k)-e(i+1,j,k+1)
695  endif
696  ! Thickness deficits can not arise simply because a layer's bottom is bounded
697  ! by the bathymetry.
698  if (e(i,j,k+1) < e(i+1,j,nz+1)) then
699  if (h1 > h2) h1 = max(e(i,j,k)-e(i+1,j,nz+1), h2)
700  elseif (e(i+1,j,k+1) < e(i,j,nz+1)) then
701  if (h2 > h1) h2 = max(e(i+1,j,k)-e(i,j,nz+1), h1)
702  endif
703  h_def_u(i,j) = h_def_u(i,j) + 0.5*(h1-h2)**2 / ((h1 + h2) + h_neglect)
704  h_norm_u(i,j) = h_norm_u(i,j) + 0.5*(h1+h2)
705  enddo ; enddo ; enddo
706  if (present(def_rat_u_2lay)) then ; do j=js,je ; do i=is-1,ie
707  def_rat_u(i,j) = g%mask2dCu(i,j) * h_def_u(i,j) / &
708  (max(hmix_min, h_norm_u(i,j)) + h_neglect)
709  def_rat_u_2lay(i,j) = g%mask2dCu(i,j) * h_def2_u(i,j) / &
710  (max(hmix_min, h_norm_u(i,j)) + h_neglect)
711  enddo ; enddo ; else ; do j=js,je ; do i=is-1,ie
712  def_rat_u(i,j) = g%mask2dCu(i,j) * h_def_u(i,j) / &
713  (max(hmix_min, h_norm_u(i,j)) + h_neglect)
714  enddo ; enddo ; endif
715 
716  ! Determine which meridional faces are problematic.
717  do j=js-1,je ; do i=is,ie
718  ! Aggregate all water below the mixed and buffer layers for the purposes of
719  ! this diagnostic.
720  h1 = e(i,j,nkmb+1)-e(i,j,nz+1) ; h2 = e(i,j+1,nkmb+1)-e(i,j+1,nz+1)
721  if (e(i,j,nz+1) < e(i,j+1,nz+1)) then
722  if (h1 > h2) h1 = max(e(i,j,nkmb+1)-e(i,j+1,nz+1), h2)
723  elseif (e(i,j+1,nz+1) < e(i,j,nz+1)) then
724  if (h2 > h1) h2 = max(e(i,j+1,nkmb+1)-e(i,j,nz+1), h1)
725  endif
726  h_def_v(i,j) = 0.5*(h1-h2)**2 / ((h1 + h2) + h_neglect)
727  h_norm_v(i,j) = 0.5*(h1+h2)
728  enddo ; enddo
729  if (present(def_rat_v_2lay)) then ; do j=js-1,je ; do i=is,ie
730  ! This is a particular metric of the aggregation into two layers.
731  h1 = e(i,j,1)-e(i,j,nkmb+1) ; h2 = e(i,j+1,1)-e(i,j+1,nkmb+1)
732  if (e(i,j,nkmb+1) < e(i,j+1,nz+1)) then
733  if (h1 > h2) h1 = max(e(i,j,1)-e(i,j+1,nz+1), h2)
734  elseif (e(i,j+1,nkmb+1) < e(i,j,nz+1)) then
735  if (h2 > h1) h2 = max(e(i,j+1,1)-e(i,j,nz+1), h1)
736  endif
737  h_def2_v(i,j) = h_def_v(i,j) + 0.5*(h1-h2)**2 / ((h1 + h2) + h_neglect)
738  enddo ; enddo ; endif
739  do k=1,nkmb ; do j=js-1,je ; do i=is,ie
740  if (present(h)) then
741  h1 = h(i,j,k) ; h2 = h(i,j+1,k)
742  else
743  h1 = e(i,j,k)-e(i,j,k+1) ; h2 = e(i,j+1,k)-e(i,j+1,k+1)
744  endif
745  ! Thickness deficits can not arise simply because a layer's bottom is bounded
746  ! by the bathymetry.
747  if (e(i,j,k+1) < e(i,j+1,nz+1)) then
748  if (h1 > h2) h1 = max(e(i,j,k)-e(i,j+1,nz+1), h2)
749  elseif (e(i,j+1,k+1) < e(i,j,nz+1)) then
750  if (h2 > h1) h2 = max(e(i,j+1,k)-e(i,j,nz+1), h1)
751  endif
752  h_def_v(i,j) = h_def_v(i,j) + 0.5*(h1-h2)**2 / ((h1 + h2) + h_neglect)
753  h_norm_v(i,j) = h_norm_v(i,j) + 0.5*(h1+h2)
754  enddo ; enddo ; enddo
755  if (present(def_rat_v_2lay)) then ; do j=js-1,je ; do i=is,ie
756  def_rat_v(i,j) = g%mask2dCv(i,j) * h_def_v(i,j) / &
757  (max(hmix_min, h_norm_v(i,j)) + h_neglect)
758  def_rat_v_2lay(i,j) = g%mask2dCv(i,j) * h_def2_v(i,j) / &
759  (max(hmix_min, h_norm_v(i,j)) + h_neglect)
760  enddo ; enddo ; else ; do j=js-1,je ; do i=is,ie
761  def_rat_v(i,j) = g%mask2dCv(i,j) * h_def_v(i,j) / &
762  (max(hmix_min, h_norm_v(i,j)) + h_neglect)
763  enddo ; enddo ; endif
764 

◆ regularize_layers()

subroutine, public mom_regularize_layers::regularize_layers ( real, dimension(szi_(g),szj_(g),szk_(g)), intent(inout)  h,
type(thermo_var_ptrs), intent(inout)  tv,
real, intent(in)  dt,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(inout)  ea,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(inout)  eb,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(regularize_layers_cs), pointer  CS 
)

This subroutine partially steps the bulk mixed layer model. The following processes are executed, in the order listed.

Parameters
[in,out]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in,out]hLayer thicknesses [H ~> m or kg m-2].
[in,out]tvA structure containing pointers to any available thermodynamic fields. Absent fields have NULL ptrs.
[in]dtTime increment [T ~> s].
[in,out]eaThe amount of fluid moved downward into a
[in,out]ebThe amount of fluid moved upward into a layer
[in]usA dimensional unit scaling type
csThe control structure returned by a previous call to regularize_layers_init.

Definition at line 72 of file MOM_regularize_layers.F90.

72  type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.
73  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
74  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
75  intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2].
76  type(thermo_var_ptrs), intent(inout) :: tv !< A structure containing pointers to any
77  !! available thermodynamic fields. Absent fields
78  !! have NULL ptrs.
79  real, intent(in) :: dt !< Time increment [T ~> s].
80  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
81  intent(inout) :: ea !< The amount of fluid moved downward into a
82  !! layer; this should be increased due to mixed
83  !! layer detrainment [H ~> m or kg m-2].
84  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
85  intent(inout) :: eb !< The amount of fluid moved upward into a layer
86  !! this should be increased due to mixed layer
87  !! entrainment [H ~> m or kg m-2].
88  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
89  type(regularize_layers_cs), pointer :: cs !< The control structure returned by a previous
90  !! call to regularize_layers_init.
91  ! Local variables
92  integer :: i, j, k, is, ie, js, je, nz
93 
94  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
95 
96  if (.not. associated(cs)) call mom_error(fatal, "MOM_regularize_layers: "//&
97  "Module must be initialized before it is used.")
98 
99  if (cs%regularize_surface_layers) then
100  call pass_var(h, g%Domain, clock=id_clock_pass)
101  call regularize_surface(h, tv, dt, ea, eb, g, gv, us, cs)
102  endif
103 

◆ regularize_layers_init()

subroutine, public mom_regularize_layers::regularize_layers_init ( type(time_type), intent(in), target  Time,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(param_file_type), intent(in)  param_file,
type(diag_ctrl), intent(inout), target  diag,
type(regularize_layers_cs), pointer  CS 
)

Initializes the regularize_layers control structure.

Parameters
[in]timeThe current model time.
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]param_fileA structure to parse for run-time parameters.
[in,out]diagA structure that is used to regulate diagnostic output.
csA pointer that is set to point to the control structure for this module.

Definition at line 769 of file MOM_regularize_layers.F90.

769  type(time_type), target, intent(in) :: time !< The current model time.
770  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
771  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
772  type(param_file_type), intent(in) :: param_file !< A structure to parse for
773  !! run-time parameters.
774  type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate
775  !! diagnostic output.
776  type(regularize_layers_cs), pointer :: cs !< A pointer that is set to point to the
777  !! control structure for this module.
778 #include "version_variable.h"
779  character(len=40) :: mdl = "MOM_regularize_layers" ! This module's name.
780  logical :: use_temperature
781  logical :: default_2018_answers
782  logical :: just_read
783  integer :: isd, ied, jsd, jed
784  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
785 
786  if (associated(cs)) then
787  call mom_error(warning, "regularize_layers_init called with an associated"// &
788  "associated control structure.")
789  return
790  else ; allocate(cs) ; endif
791 
792  cs%diag => diag
793  cs%Time => time
794 
795 ! Set default, read and log parameters
796  call get_param(param_file, mdl, "REGULARIZE_SURFACE_LAYERS", cs%regularize_surface_layers, &
797  default=.false., do_not_log=.true.)
798  call log_version(param_file, mdl, version, "", all_default=.not.cs%regularize_surface_layers)
799  call get_param(param_file, mdl, "REGULARIZE_SURFACE_LAYERS", cs%regularize_surface_layers, &
800  "If defined, vertically restructure the near-surface "//&
801  "layers when they have too much lateral variations to "//&
802  "allow for sensible lateral barotropic transports.", &
803  default=.false.)
804  just_read = .not.cs%regularize_surface_layers
805  if (cs%regularize_surface_layers) then
806  call get_param(param_file, mdl, "REGULARIZE_SURFACE_DETRAIN", cs%reg_sfc_detrain, &
807  "If true, allow the buffer layers to detrain into the "//&
808  "interior as a part of the restructuring when "//&
809  "REGULARIZE_SURFACE_LAYERS is true.", default=.true., do_not_log=just_read)
810  call get_param(param_file, mdl, "REG_SFC_DENSE_MATCH_TOLERANCE", cs%density_match_tol, &
811  "A relative tolerance for how well the densities must match with the target "//&
812  "densities during detrainment when regularizing the near-surface layers. The "//&
813  "default of 0.6 gives 20% overlaps in density", &
814  units="nondim", default=0.6, do_not_log=just_read)
815  call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, &
816  "This sets the default value for the various _2018_ANSWERS parameters.", &
817  default=.false., do_not_log=just_read)
818  call get_param(param_file, mdl, "REGULARIZE_LAYERS_2018_ANSWERS", cs%answers_2018, &
819  "If true, use the order of arithmetic and expressions that recover the answers "//&
820  "from the end of 2018. Otherwise, use updated and more robust forms of the "//&
821  "same expressions.", default=default_2018_answers, do_not_log=just_read)
822  endif
823 
824  call get_param(param_file, mdl, "HMIX_MIN", cs%Hmix_min, &
825  "The minimum mixed layer depth if the mixed layer depth is determined "//&
826  "dynamically.", units="m", default=0.0, scale=gv%m_to_H, do_not_log=just_read)
827  call get_param(param_file, mdl, "REG_SFC_DEFICIT_TOLERANCE", cs%h_def_tol1, &
828  "The value of the relative thickness deficit at which "//&
829  "to start modifying the layer structure when "//&
830  "REGULARIZE_SURFACE_LAYERS is true.", units="nondim", &
831  default=0.5, do_not_log=just_read)
832  cs%h_def_tol2 = 0.2 + 0.8*cs%h_def_tol1
833  cs%h_def_tol3 = 0.3 + 0.7*cs%h_def_tol1
834  cs%h_def_tol4 = 0.5 + 0.5*cs%h_def_tol1
835 
836  call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false.)
837 ! if (.not. CS%debug) &
838 ! call get_param(param_file, mdl, "DEBUG_CONSERVATION", CS%debug, &
839 ! "If true, monitor conservation and extrema.", default=.false., do_not_log=just_read)
840 
841  call get_param(param_file, mdl, "ALLOW_CLOCKS_IN_OMP_LOOPS", cs%allow_clocks_in_omp_loops, &
842  "If true, clocks can be called from inside loops that can "//&
843  "be threaded. To run with multiple threads, set to False.", &
844  default=.true., do_not_log=just_read)
845 
846  if (.not.cs%regularize_surface_layers) return
847 
848  cs%id_def_rat = register_diag_field('ocean_model', 'deficit_ratio', diag%axesT1, &
849  time, 'Max face thickness deficit ratio', 'nondim')
850 
851  if (cs%allow_clocks_in_omp_loops) then
852  id_clock_eos = cpu_clock_id('(Ocean regularize_layers EOS)', grain=clock_routine)
853  endif
854  id_clock_pass = cpu_clock_id('(Ocean regularize_layers halo updates)', grain=clock_routine)
855 

◆ regularize_surface()

subroutine mom_regularize_layers::regularize_surface ( real, dimension(szi_(g),szj_(g),szk_(g)), intent(inout)  h,
type(thermo_var_ptrs), intent(inout)  tv,
real, intent(in)  dt,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(inout)  ea,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(inout)  eb,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(regularize_layers_cs), pointer  CS 
)
private

This subroutine ensures that there is a degree of horizontal smoothness in the depths of the near-surface interfaces.

Parameters
[in,out]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in,out]hLayer thicknesses [H ~> m or kg m-2].
[in,out]tvA structure containing pointers to any available thermodynamic fields. Absent fields have NULL ptrs.
[in]dtTime increment [T ~> s].
[in,out]eaThe amount of fluid moved downward into a
[in,out]ebThe amount of fluid moved upward into a layer
[in]usA dimensional unit scaling type
csThe control structure returned by a previous call to regularize_layers_init.

Definition at line 109 of file MOM_regularize_layers.F90.

109  type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.
110  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
111  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
112  intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2].
113  type(thermo_var_ptrs), intent(inout) :: tv !< A structure containing pointers to any
114  !! available thermodynamic fields. Absent fields
115  !! have NULL ptrs.
116  real, intent(in) :: dt !< Time increment [T ~> s].
117  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
118  intent(inout) :: ea !< The amount of fluid moved downward into a
119  !! layer; this should be increased due to mixed
120  !! layer detrainment [H ~> m or kg m-2].
121  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
122  intent(inout) :: eb !< The amount of fluid moved upward into a layer
123  !! this should be increased due to mixed layer
124  !! entrainment [H ~> m or kg m-2].
125  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
126  type(regularize_layers_cs), pointer :: cs !< The control structure returned by a previous
127  !! call to regularize_layers_init.
128  ! Local variables
129  real, dimension(SZIB_(G),SZJ_(G)) :: &
130  def_rat_u ! The ratio of the thickness deficit to the minimum depth [nondim].
131  real, dimension(SZI_(G),SZJB_(G)) :: &
132  def_rat_v ! The ratio of the thickness deficit to the minimum depth [nondim].
133  real, dimension(SZI_(G),SZJ_(G)) :: &
134  def_rat_h ! The ratio of the thickness deficit to the minimum depth [nondim].
135  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: &
136  e ! The interface depths [H ~> m or kg m-2], positive upward.
137 
138  real, dimension(SZI_(G),SZK_(G)+1) :: &
139  e_filt, e_2d ! The interface depths [H ~> m or kg m-2], positive upward.
140  real, dimension(SZI_(G),SZK_(G)) :: &
141  h_2d, & ! A 2-d version of h [H ~> m or kg m-2].
142  t_2d, & ! A 2-d version of tv%T [degC].
143  s_2d, & ! A 2-d version of tv%S [ppt].
144  rcv, & ! A 2-d version of the coordinate density [R ~> kg m-3].
145  h_2d_init, & ! The initial value of h_2d [H ~> m or kg m-2].
146  t_2d_init, & ! THe initial value of T_2d [degC].
147  s_2d_init, & ! The initial value of S_2d [ppt].
148  d_eb, & ! The downward increase across a layer in the entrainment from
149  ! below [H ~> m or kg m-2]. The sign convention is that positive values of
150  ! d_eb correspond to a gain in mass by a layer by upward motion.
151  d_ea ! The upward increase across a layer in the entrainment from
152  ! above [H ~> m or kg m-2]. The sign convention is that positive values of
153  ! d_ea mean a net gain in mass by a layer from downward motion.
154  real, dimension(SZI_(G)) :: &
155  p_ref_cv, & ! Reference pressure for the potential density which defines
156  ! the coordinate variable, set to P_Ref [R L2 T-2 ~> Pa].
157  rcv_tol, & ! A tolerence, relative to the target density differences
158  ! between layers, for detraining into the interior [nondim].
159  h_add_tgt, h_add_tot, &
160  h_tot1, th_tot1, sh_tot1, &
161  h_tot3, th_tot3, sh_tot3, &
162  h_tot2, th_tot2, sh_tot2
163  real, dimension(SZK_(G)) :: &
164  h_prev_1d ! The previous thicknesses [H ~> m or kg m-2].
165  real :: i_dtol ! The inverse of the tolerance changes [nondim].
166  real :: i_dtol34 ! The inverse of the tolerance changes [nondim].
167  real :: h1, h2 ! Temporary thicknesses [H ~> m or kg m-2].
168  real :: e_e, e_w, e_n, e_s ! Temporary interface heights [H ~> m or kg m-2].
169  real :: wt ! The weight of the filted interfaces in setting the targets [nondim].
170  real :: scale ! A scaling factor [nondim].
171  real :: h_neglect ! A thickness that is so small it is usually lost
172  ! in roundoff and can be neglected [H ~> m or kg m-2].
173  real, dimension(SZK_(G)+1) :: &
174  int_flux, int_tflux, int_sflux, int_rflux
175  real :: h_add
176  real :: h_det_tot
177  real :: max_def_rat
178  real :: rcv_min_det ! The lightest (min) and densest (max) coordinate density
179  real :: rcv_max_det ! that can detrain into a layer [R ~> kg m-3].
180 
181  real :: int_top, int_bot
182  real :: h_predicted
183  real :: h_prev
184  real :: h_deficit
185 
186  logical :: cols_left, ent_any, more_ent_i(szi_(g)), ent_i(szi_(g))
187  logical :: det_any, det_i(szi_(g))
188  logical :: do_j(szj_(g)), do_i(szi_(g)), find_i(szi_(g))
189  logical :: debug = .false.
190  logical :: fatal_error
191  character(len=256) :: mesg ! Message for error messages.
192  integer, dimension(2) :: eosdom ! The i-computational domain for the equation of state
193  integer :: i, j, k, is, ie, js, je, nz, nkmb, nkml, k1, k2, k3, ks, nz_filt, kmax_d_ea
194 
195  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
196 
197  if (.not. associated(cs)) call mom_error(fatal, "MOM_regularize_layers: "//&
198  "Module must be initialized before it is used.")
199 
200  if (gv%nkml<1) return
201  nkmb = gv%nk_rho_varies ; nkml = gv%nkml
202  if (.not.associated(tv%eqn_of_state)) call mom_error(fatal, &
203  "MOM_regularize_layers: This module now requires the use of temperature and "//&
204  "an equation of state.")
205 
206  h_neglect = gv%H_subroundoff
207  debug = (debug .or. cs%debug)
208 
209  i_dtol = 1.0 / max(cs%h_def_tol2 - cs%h_def_tol1, 1e-40)
210  i_dtol34 = 1.0 / max(cs%h_def_tol4 - cs%h_def_tol3, 1e-40)
211 
212  p_ref_cv(:) = tv%P_Ref
213  eosdom(:) = eos_domain(g%HI)
214 
215  do j=js-1,je+1 ; do i=is-1,ie+1
216  e(i,j,1) = 0.0
217  enddo ; enddo
218  do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
219  e(i,j,k+1) = e(i,j,k) - h(i,j,k)
220  enddo ; enddo ; enddo
221 
222  call find_deficit_ratios(e, def_rat_u, def_rat_v, g, gv, cs, h=h)
223 
224  ! Determine which columns are problematic
225  do j=js,je ; do_j(j) = .false. ; enddo
226  do j=js,je ; do i=is,ie
227  def_rat_h(i,j) = max(def_rat_u(i-1,j), def_rat_u(i,j), &
228  def_rat_v(i,j-1), def_rat_v(i,j))
229  if (def_rat_h(i,j) > cs%h_def_tol1) do_j(j) = .true.
230  enddo ; enddo
231 
232  ! Now restructure the layers.
233  !$OMP parallel do default(private) shared(is,ie,js,je,nz,do_j,def_rat_h,CS,nkmb,G,GV,US, &
234  !$OMP e,I_dtol,h,tv,debug,h_neglect,p_ref_cv,ea, &
235  !$OMP eb,id_clock_EOS,nkml,EOSdom)
236  do j=js,je ; if (do_j(j)) then
237 
238 ! call cpu_clock_begin(id_clock_EOS)
239 ! call calculate_density_derivs(T(:,1), S(:,1), p_ref_cv, dRcv_dT, dRcv_dS, tv%eqn_of_state, EOSdom)
240 ! call cpu_clock_end(id_clock_EOS)
241 
242  do k=1,nz ; do i=is,ie ; d_ea(i,k) = 0.0 ; d_eb(i,k) = 0.0 ; enddo ; enddo
243  kmax_d_ea = 0
244 
245  max_def_rat = 0.0
246  do i=is,ie
247  do_i(i) = def_rat_h(i,j) > cs%h_def_tol1
248  if (def_rat_h(i,j) > max_def_rat) max_def_rat = def_rat_h(i,j)
249  enddo
250  nz_filt = nkmb+1 ; if (max_def_rat > cs%h_def_tol3) nz_filt = nz+1
251 
252  ! Find a 2-D 1-2-1 filtered version of e to target. Area weights are
253  ! deliberately omitted here. This is slightly more complicated than a
254  ! simple filter so that the effects of topography are eliminated.
255  do k=1,nz_filt ; do i=is,ie ; if (do_i(i)) then
256  if (g%mask2dCu(i,j) <= 0.0) then ; e_e = e(i,j,k) ; else
257  e_e = max(e(i+1,j,k) + min(e(i,j,k) - e(i+1,j,nz+1), 0.0), &
258  e(i,j,nz+1) + (nz+1-k)*gv%Angstrom_H)
259 
260  endif
261  if (g%mask2dCu(i-1,j) <= 0.0) then ; e_w = e(i,j,k) ; else
262  e_w = max(e(i-1,j,k) + min(e(i,j,k) - e(i-1,j,nz+1), 0.0), &
263  e(i,j,nz+1) + (nz+1-k)*gv%Angstrom_H)
264  endif
265  if (g%mask2dCv(i,j) <= 0.0) then ; e_n = e(i,j,k) ; else
266  e_n = max(e(i,j+1,k) + min(e(i,j,k) - e(i,j+1,nz+1), 0.0), &
267  e(i,j,nz+1) + (nz+1-k)*gv%Angstrom_H)
268  endif
269  if (g%mask2dCv(i,j-1) <= 0.0) then ; e_s = e(i,j,k) ; else
270  e_s = max(e(i,j-1,k) + min(e(i,j,k) - e(i,j-1,nz+1), 0.0), &
271  e(i,j,nz+1) + (nz+1-k)*gv%Angstrom_H)
272  endif
273 
274  wt = max(0.0, min(1.0, i_dtol*(def_rat_h(i,j)-cs%h_def_tol1)))
275 
276  e_filt(i,k) = (1.0 - 0.5*wt) * e(i,j,k) + &
277  wt * 0.125 * ((e_e + e_w) + (e_n + e_s))
278  e_2d(i,k) = e(i,j,k)
279  endif ; enddo ; enddo
280  do k=1,nz ; do i=is,ie
281  h_2d(i,k) = h(i,j,k)
282  t_2d(i,k) = tv%T(i,j,k) ; s_2d(i,k) = tv%S(i,j,k)
283  enddo ; enddo
284 
285  if (debug) then
286  do k=1,nz ; do i=is,ie ; if (do_i(i)) then
287  h_2d_init(i,k) = h(i,j,k)
288  t_2d_init(i,k) = tv%T(i,j,k) ; s_2d_init(i,k) = tv%S(i,j,k)
289  endif ; enddo ; enddo
290  endif
291 
292  ! First, try to entrain from the interior.
293  ent_any = .false.
294  do i=is,ie
295  more_ent_i(i) = .false. ; ent_i(i) = .false.
296  h_add_tgt(i) = 0.0 ; h_add_tot(i) = 0.0
297  if (do_i(i) .and. (e_2d(i,nkmb+1) > e_filt(i,nkmb+1))) then
298  more_ent_i(i) = .true. ; ent_i(i) = .true. ; ent_any = .true.
299  h_add_tgt(i) = e_2d(i,nkmb+1) - e_filt(i,nkmb+1)
300  endif
301  enddo
302 
303  if (ent_any) then
304  do k=nkmb+1,nz
305  cols_left = .false.
306  do i=is,ie ; if (more_ent_i(i)) then
307  if (h_2d(i,k) - gv%Angstrom_H > h_neglect) then
308  if (e_2d(i,nkmb+1)-e_filt(i,nkmb+1) > h_2d(i,k) - gv%Angstrom_H) then
309  h_add = h_2d(i,k) - gv%Angstrom_H
310  h_2d(i,k) = gv%Angstrom_H
311  e_2d(i,nkmb+1) = e_2d(i,nkmb+1) - h_add
312  else
313  h_add = e_2d(i,nkmb+1) - e_filt(i,nkmb+1)
314  h_2d(i,k) = h_2d(i,k) - h_add
315  if (cs%answers_2018) then
316  e_2d(i,nkmb+1) = e_2d(i,nkmb+1) - h_add
317  else
318  e_2d(i,nkmb+1) = e_filt(i,nkmb+1)
319  endif
320  endif
321  d_eb(i,k-1) = d_eb(i,k-1) + h_add
322  h_add_tot(i) = h_add_tot(i) + h_add
323  h_prev = h_2d(i,nkmb)
324  h_2d(i,nkmb) = h_2d(i,nkmb) + h_add
325 
326  t_2d(i,nkmb) = (h_prev*t_2d(i,nkmb) + h_add*t_2d(i,k)) / h_2d(i,nkmb)
327  s_2d(i,nkmb) = (h_prev*s_2d(i,nkmb) + h_add*s_2d(i,k)) / h_2d(i,nkmb)
328 
329  if ((e_2d(i,nkmb+1) <= e_filt(i,nkmb+1)) .or. &
330  (h_add_tot(i) > 0.6*h_add_tgt(i))) then !### 0.6 is adjustable?.
331  more_ent_i(i) = .false.
332  else
333  cols_left = .true.
334  endif
335  else
336  cols_left = .true.
337  endif
338  endif ; enddo
339  if (.not.cols_left) exit
340  enddo
341 
342  ks = min(k-1,nz-1)
343  do k=ks,nkmb,-1 ; do i=is,ie ; if (ent_i(i)) then
344  d_eb(i,k) = d_eb(i,k) + d_eb(i,k+1)
345  endif ; enddo ; enddo
346  endif ! ent_any
347 
348  ! This is where code to detrain to the interior will go.
349  ! The buffer layers can only detrain water into layers when the buffer
350  ! layer potential density is between (c*Rlay(k-1) + (1-c)*Rlay(k)) and
351  ! (c*Rlay(k+1) + (1-c)*Rlay(k)), where 0.5 <= c < 1.0.
352  ! Do not detrain if the 2-layer deficit ratio is not significant.
353  ! Detrainment must be able to come from all mixed and buffer layers.
354  ! All water is moved out of the buffer layers below before moving from
355  ! a shallower layer (characteristics do not cross).
356  det_any = .false.
357  if ((max_def_rat > cs%h_def_tol3) .and. (cs%reg_sfc_detrain)) then
358  do i=is,ie
359  det_i(i) = .false. ; rcv_tol(i) = 0.0
360  if (do_i(i) .and. (e_2d(i,nkmb+1) < e_filt(i,nkmb+1)) .and. &
361  (def_rat_h(i,j) > cs%h_def_tol3)) then
362  det_i(i) = .true. ; det_any = .true.
363  ! The CS%density_match_tol default value of 0.6 gives 20% overlap in acceptable densities.
364  rcv_tol(i) = cs%density_match_tol * min((def_rat_h(i,j) - cs%h_def_tol3), 1.0)
365  endif
366  enddo
367  endif
368  if (det_any) then
369  call cpu_clock_begin(id_clock_eos)
370  do k=1,nkmb
371  call calculate_density(t_2d(:,k), s_2d(:,k), p_ref_cv, rcv(:,k), tv%eqn_of_state, eosdom)
372  enddo
373  call cpu_clock_end(id_clock_eos)
374 
375  do i=is,ie ; if (det_i(i)) then
376  k1 = nkmb ; k2 = nz
377  h_det_tot = 0.0
378  do ! This loop is terminated by exits.
379  if (k1 <= 1) exit
380  if (k2 <= nkmb) exit
381  rcv_min_det = (gv%Rlay(k2) + rcv_tol(i)*(gv%Rlay(k2-1)-gv%Rlay(k2)))
382  if (k2 < nz) then
383  rcv_max_det = (gv%Rlay(k2) + rcv_tol(i)*(gv%Rlay(k2+1)-gv%Rlay(k2)))
384  else
385  rcv_max_det = (gv%Rlay(nz) + rcv_tol(i)*(gv%Rlay(nz)-gv%Rlay(nz-1)))
386  endif
387  if (rcv(i,k1) > rcv_max_det) &
388  exit ! All shallower interior layers are too light for detrainment.
389 
390  h_deficit = (e_filt(i,k2)-e_filt(i,k2+1)) - h_2d(i,k2)
391  if ((e_filt(i,k2) > e_2d(i,k1+1)) .and. (h_deficit > 0.0) .and. &
392  (rcv(i,k1) < rcv_max_det) .and. (rcv(i,k1) > rcv_min_det)) then
393  ! Detrainment will occur.
394  h_add = min(e_filt(i,k2) - e_2d(i,k2), h_deficit )
395  if (h_add < h_2d(i,k1)) then
396  ! Only part of layer k1 detrains.
397  if (h_add > 0.0) then
398  h_prev = h_2d(i,k2)
399  h_2d(i,k2) = h_2d(i,k2) + h_add
400  e_2d(i,k2) = e_2d(i,k2+1) + h_2d(i,k2)
401  d_ea(i,k2) = d_ea(i,k2) + h_add
402  kmax_d_ea = max(kmax_d_ea, k2)
403  ! This is upwind. It should perhaps be higher order...
404  t_2d(i,k2) = (h_prev*t_2d(i,k2) + h_add*t_2d(i,k1)) / h_2d(i,k2)
405  s_2d(i,k2) = (h_prev*s_2d(i,k2) + h_add*s_2d(i,k1)) / h_2d(i,k2)
406  h_det_tot = h_det_tot + h_add
407 
408  h_2d(i,k1) = h_2d(i,k1) - h_add
409  do k3=k1,nkmb ; e_2d(i,k3+1) = e_2d(i,k3) - h_2d(i,k3) ; enddo
410  do k3=k1+1,nkmb ; d_ea(i,k3) = d_ea(i,k3) + h_add ; enddo
411  else
412  if (h_add < 0.0) &
413  call mom_error(fatal, "h_add is negative. Some logic is wrong.")
414  h_add = 0.0 ! This usually should not happen...
415  endif
416 
417  ! Move up to the next target layer.
418  k2 = k2-1
419  if (k2>nkmb+1) e_2d(i,k2) = e_2d(i,k2) + h_det_tot
420  else
421  h_add = h_2d(i,k1)
422  h_prev = h_2d(i,k2)
423  h_2d(i,k2) = h_2d(i,k2) + h_add
424  e_2d(i,k2) = e_2d(i,k2+1) + h_2d(i,k2)
425  d_ea(i,k2) = d_ea(i,k2) + h_add
426  kmax_d_ea = max(kmax_d_ea, k2)
427  t_2d(i,k2) = (h_prev*t_2d(i,k2) + h_add*t_2d(i,k1)) / h_2d(i,k2)
428  s_2d(i,k2) = (h_prev*s_2d(i,k2) + h_add*s_2d(i,k1)) / h_2d(i,k2)
429  h_det_tot = h_det_tot + h_add
430 
431  h_2d(i,k1) = 0.0
432  do k3=k1,nkmb ; e_2d(i,k3+1) = e_2d(i,k3) - h_2d(i,k3) ; enddo
433  do k3=k1+1,nkmb ; d_ea(i,k3) = d_ea(i,k3) + h_add ; enddo
434 
435  ! Move up to the next source layer.
436  k1 = k1-1
437  endif
438 
439  else
440  ! Move up to the next target layer.
441  k2 = k2-1
442  if (k2>nkmb+1) e_2d(i,k2) = e_2d(i,k2) + h_det_tot
443  endif
444 
445  enddo ! exit terminated loop.
446  endif ; enddo
447  do k=kmax_d_ea-1,nkmb+1,-1 ; do i=is,ie ; if (det_i(i)) then
448  d_ea(i,k) = d_ea(i,k) + d_ea(i,k+1)
449  endif ; enddo ; enddo
450  endif ! Detrainment to the interior.
451  if (debug) then
452  do i=is,ie ; h_tot3(i) = 0.0 ; th_tot3(i) = 0.0 ; sh_tot3(i) = 0.0 ; enddo
453  do k=1,nz ; do i=is,ie ; if (do_i(i)) then
454  h_tot3(i) = h_tot3(i) + h_2d(i,k)
455  th_tot3(i) = th_tot3(i) + h_2d(i,k) * t_2d(i,k)
456  sh_tot3(i) = sh_tot3(i) + h_2d(i,k) * s_2d(i,k)
457  endif ; enddo ; enddo
458  endif
459 
460  do i=is,ie ; if (do_i(i)) then
461  ! Rescale the interface targets so the depth at the bottom of the deepest
462  ! buffer layer matches.
463  scale = e_2d(i,nkmb+1) / e_filt(i,nkmb+1)
464  do k=2,nkmb+1 ; e_filt(i,k) = e_filt(i,k) * scale ; enddo
465 
466  ! Ensure that layer 1 only has water from layers 1 to nkml and rescale
467  ! the remaining layer thicknesses if necessary.
468  if (e_filt(i,2) < e_2d(i,nkml)) then
469  scale = (e_2d(i,nkml) - e_filt(i,nkmb+1)) / &
470  ((e_filt(i,2) - e_filt(i,nkmb+1)) + h_neglect)
471  do k=3,nkmb
472  e_filt(i,k) = e_filt(i,nkmb+1) + scale * (e_filt(i,k) - e_filt(i,nkmb+1))
473  enddo
474  e_filt(i,2) = e_2d(i,nkml)
475  endif
476 
477  ! Map the water back into the layers. There are not mixed or buffer layers that are exceedingly
478  ! small compared to the others, so the code here is less prone to roundoff than elsewhere in MOM6.
479  k1 = 1 ; k2 = 1
480  int_top = 0.0
481  do k=1,nkmb+1
482  int_flux(k) = 0.0 ; int_rflux(k) = 0.0
483  int_tflux(k) = 0.0 ; int_sflux(k) = 0.0
484  enddo
485  do k=1,2*nkmb
486  int_bot = max(e_2d(i,k1+1),e_filt(i,k2+1))
487  h_add = int_top - int_bot
488 
489  if (k2 > k1) then
490  do k3=k1+1,k2
491  d_ea(i,k3) = d_ea(i,k3) + h_add
492  int_flux(k3) = int_flux(k3) + h_add
493  int_tflux(k3) = int_tflux(k3) + h_add*t_2d(i,k1)
494  int_sflux(k3) = int_sflux(k3) + h_add*s_2d(i,k1)
495  enddo
496  elseif (k1 > k2) then
497  do k3=k2,k1-1
498  d_eb(i,k3) = d_eb(i,k3) + h_add
499  int_flux(k3+1) = int_flux(k3+1) - h_add
500  int_tflux(k3+1) = int_tflux(k3+1) - h_add*t_2d(i,k1)
501  int_sflux(k3+1) = int_sflux(k3+1) - h_add*s_2d(i,k1)
502  enddo
503  endif
504 
505  if (int_bot <= e_filt(i,k2+1)) then
506  ! Increment the target layer.
507  k2 = k2 + 1
508  elseif (int_bot <= e_2d(i,k1+1)) then
509  ! Increment the source layer.
510  k1 = k1 + 1
511  else
512  call mom_error(fatal, &
513  "Regularize_surface: Could not increment target or source.")
514  endif
515  if ((k1 > nkmb) .or. (k2 > nkmb)) exit
516  int_top = int_bot
517  enddo
518  if (k2 < nkmb) &
519  call mom_error(fatal, "Regularize_surface: Did not assign fluid to layer nkmb.")
520 
521  ! Note that movement of water across the base of the bottommost buffer
522  ! layer has already been dealt with separately.
523  do k=1,nkmb ; h_prev_1d(k) = h_2d(i,k) ; enddo
524  h_2d(i,1) = h_2d(i,1) - int_flux(2)
525  do k=2,nkmb-1
526  h_2d(i,k) = h_2d(i,k) + (int_flux(k) - int_flux(k+1))
527  enddo
528  ! Note that movement of water across the base of the bottommost buffer
529  ! layer has already been dealt with separately.
530  h_2d(i,nkmb) = h_2d(i,nkmb) + int_flux(nkmb)
531 
532  t_2d(i,1) = (t_2d(i,1)*h_prev_1d(1) - int_tflux(2)) / h_2d(i,1)
533  s_2d(i,1) = (s_2d(i,1)*h_prev_1d(1) - int_sflux(2)) / h_2d(i,1)
534  do k=2,nkmb-1
535  t_2d(i,k) = (t_2d(i,k)*h_prev_1d(k) + (int_tflux(k) - int_tflux(k+1))) / h_2d(i,k)
536  s_2d(i,k) = (s_2d(i,k)*h_prev_1d(k) + (int_sflux(k) - int_sflux(k+1))) / h_2d(i,k)
537  enddo
538  t_2d(i,nkmb) = (t_2d(i,nkmb)*h_prev_1d(nkmb) + int_tflux(nkmb) ) / h_2d(i,nkmb)
539  s_2d(i,nkmb) = (s_2d(i,nkmb)*h_prev_1d(nkmb) + int_sflux(nkmb) ) / h_2d(i,nkmb)
540 
541  endif ; enddo ! i-loop
542 
543  ! Copy the interior thicknesses and other fields back to the 3-d arrays.
544  do k=1,nz ; do i=is,ie ; if (do_i(i)) then
545  h(i,j,k) = h_2d(i,k)
546  tv%T(i,j,k) = t_2d(i,k) ; tv%S(i,j,k) = s_2d(i,k)
547  ea(i,j,k) = ea(i,j,k) + d_ea(i,k)
548  eb(i,j,k) = eb(i,j,k) + d_eb(i,k)
549  endif ; enddo ; enddo
550 
551  if (debug) then
552  do i=is,ie ; h_tot1(i) = 0.0 ; th_tot1(i) = 0.0 ; sh_tot1(i) = 0.0 ; enddo
553  do i=is,ie ; h_tot2(i) = 0.0 ; th_tot2(i) = 0.0 ; sh_tot2(i) = 0.0 ; enddo
554 
555  do k=1,nz ; do i=is,ie ; if (do_i(i)) then
556  h_tot1(i) = h_tot1(i) + h_2d_init(i,k)
557  h_tot2(i) = h_tot2(i) + h(i,j,k)
558 
559  th_tot1(i) = th_tot1(i) + h_2d_init(i,k) * t_2d_init(i,k)
560  th_tot2(i) = th_tot2(i) + h(i,j,k) * tv%T(i,j,k)
561  sh_tot1(i) = sh_tot1(i) + h_2d_init(i,k) * s_2d_init(i,k)
562  sh_tot2(i) = sh_tot2(i) + h(i,j,k) * tv%S(i,j,k)
563  if (h(i,j,k) < 0.0) &
564  call mom_error(fatal,"regularize_surface: Negative thicknesses.")
565  if (k==1) then ; h_predicted = h_2d_init(i,k) + (d_eb(i,k) - d_ea(i,k+1))
566  elseif (k==nz) then ; h_predicted = h_2d_init(i,k) + (d_ea(i,k) - d_eb(i,k-1))
567  else
568  h_predicted = h_2d_init(i,k) + ((d_ea(i,k) - d_eb(i,k-1)) + &
569  (d_eb(i,k) - d_ea(i,k+1)))
570  endif
571  if (abs(h(i,j,k) - h_predicted) > max(1e-9*abs(h_predicted),gv%Angstrom_H)) &
572  call mom_error(fatal, "regularize_surface: d_ea mismatch.")
573  endif ; enddo ; enddo
574  do i=is,ie ; if (do_i(i)) then
575  fatal_error = .false.
576  if (abs(h_tot1(i) - h_tot2(i)) > 1e-12*h_tot1(i)) then
577  write(mesg,'(ES11.4," became ",ES11.4," diff ",ES11.4)') &
578  h_tot1(i), h_tot2(i), (h_tot1(i) - h_tot2(i))
579  call mom_error(warning, "regularize_surface: Mass non-conservation."//&
580  trim(mesg), .true.)
581  fatal_error = .true.
582  endif
583  if (abs(th_tot1(i) - th_tot2(i)) > 1e-12*(th_tot1(i)+10.0*h_tot1(i))) then
584  write(mesg,'(ES11.4," became ",ES11.4," diff ",ES11.4," int diff ",ES11.4)') &
585  th_tot1(i), th_tot2(i), (th_tot1(i) - th_tot2(i)), (th_tot1(i) - th_tot3(i))
586  call mom_error(warning, "regularize_surface: Heat non-conservation."//&
587  trim(mesg), .true.)
588  fatal_error = .true.
589  endif
590  if (abs(sh_tot1(i) - sh_tot2(i)) > 1e-12*(sh_tot1(i)+10.0*h_tot1(i))) then
591  write(mesg,'(ES11.4," became ",ES11.4," diff ",ES11.4," int diff ",ES11.4)') &
592  sh_tot1(i), sh_tot2(i), (sh_tot1(i) - sh_tot2(i)), (sh_tot1(i) - sh_tot3(i))
593  call mom_error(warning, "regularize_surface: Salinity non-conservation."//&
594  trim(mesg), .true.)
595  fatal_error = .true.
596  endif
597  if (fatal_error) then
598  write(mesg,'("Error at lat/lon ",2(ES11.4))') g%geoLatT(i,j), g%geoLonT(i,j)
599  call mom_error(fatal, "regularize_surface: Terminating with fatal error. "//&
600  trim(mesg))
601  endif
602  endif ; enddo
603  endif
604 
605  endif ; enddo ! j-loop.
606 
607  if (cs%id_def_rat > 0) call post_data(cs%id_def_rat, def_rat_h, cs%diag)
608 

Variable Documentation

◆ id_clock_eos

integer mom_regularize_layers::id_clock_eos
private

Clock IDs.

Todo:
Should these be global?

Definition at line 64 of file MOM_regularize_layers.F90.

◆ id_clock_pass

integer mom_regularize_layers::id_clock_pass

Clock IDs.

Todo:
Should these be global?

Definition at line 64 of file MOM_regularize_layers.F90.

64 integer :: id_clock_pass, id_clock_eos