MOM6
Data Types | Functions/Subroutines | Variables
mom_regularize_layers Module Reference

Data Types

type  regularize_layers_cs
 

Functions/Subroutines

subroutine, public regularize_layers (h, tv, dt, ea, eb, G, GV, 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, 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, param_file, diag, CS)
 

Variables

integer id_clock_pass
 
integer id_clock_eos
 

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, in 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, in H (usually m or kg

Definition at line 809 of file MOM_regularize_layers.F90.

Referenced by regularize_surface().

809  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
810  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
811  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), &
812  intent(in) :: e !< Interface depths, in m or kg m-2.
813  real, dimension(SZIB_(G),SZJ_(G)), &
814  intent(out) :: def_rat_u !< The thickness deficit ratio at u points,
815  !! nondim.
816  real, dimension(SZI_(G),SZJB_(G)), &
817  intent(out) :: def_rat_v !< The thickness deficit ratio at v points,
818  !! nondim.
819  type(regularize_layers_cs), pointer :: cs !< The control structure returned by a
820  !! previous call to regularize_layers_init.
821  real, dimension(SZIB_(G),SZJ_(G)), &
822  optional, intent(out) :: def_rat_u_2lay !< The thickness deficit ratio at u
823  !! points when the mixed and buffer layers
824  !! are aggregated into 1 layer, nondim.
825  real, dimension(SZI_(G),SZJB_(G)), &
826  optional, intent(out) :: def_rat_v_2lay !< The thickness deficit ratio at v
827  !! pointswhen the mixed and buffer layers
828  !! are aggregated into 1 layer, nondim.
829  integer, optional, intent(in) :: halo !< An extra-wide halo size, 0 by default.
830  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
831  optional, intent(in) :: h !< Layer thicknesses, in H (usually m or kg
832  !! m-2); if h is not present, vertical
833  !! differences in interface heights are used
834  !! instead.
835 
836 ! This subroutine determines the amount by which the harmonic mean
837 ! thickness at velocity points differ from the arithmetic means, relative to
838 ! the the arithmetic means, after eliminating thickness variations that are
839 ! solely due to topography and aggregating all interior layers into one.
840 
841 ! Arguments: e - Interface depths, in m or kg m-2.
842 ! (out) def_rat_u - The thickness deficit ratio at u points, nondim.
843 ! (out) def_rat_v - The thickness deficit ratio at v points, nondim.
844 ! (in) G - The ocean's grid structure.
845 ! (in) GV - The ocean's vertical grid structure.
846 ! (in) CS - The control structure returned by a previous call to
847 ! regularize_layers_init.
848 ! (out,opt) def_rat_u_2lay - The thickness deficit ratio at u points when the
849 ! mixed and buffer layers are aggregated into 1 layer, nondim.
850 ! (out,opt) def_rat_v_2lay - The thickness deficit ratio at v pointswhen the
851 ! mixed and buffer layers are aggregated into 1 layer, nondim.
852 ! (in,opt) halo - An extra-wide halo size, 0 by default.
853 ! (in,opt) h - The layer thicknesse; if not present take vertical differences of e.
854  real, dimension(SZIB_(G),SZJ_(G)) :: &
855  h_def_u, & ! The vertically summed thickness deficits at u-points, in H.
856  h_norm_u, & ! The vertically summed arithmetic mean thickness by which
857  ! h_def_u is normalized, in H.
858  h_def2_u
859  real, dimension(SZI_(G),SZJB_(G)) :: &
860  h_def_v, & ! The vertically summed thickness deficits at v-points, in H.
861  h_norm_v, & ! The vertically summed arithmetic mean thickness by which
862  ! h_def_v is normalized, in H.
863  h_def2_v
864  real :: h_neglect ! A thickness that is so small it is usually lost
865  ! in roundoff and can be neglected, in H.
866  real :: h1, h2 ! Temporary thicknesses, in H.
867  integer :: i, j, k, is, ie, js, je, nz, nkmb
868 
869  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
870  if (present(halo)) then
871  is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
872  endif
873  nkmb = gv%nk_rho_varies
874  h_neglect = gv%H_subroundoff
875 
876  ! Determine which zonal faces are problematic.
877  do j=js,je ; do i=is-1,ie
878  ! Aggregate all water below the mixed and buffer layers for the purposes of
879  ! this diagnostic.
880  h1 = e(i,j,nkmb+1)-e(i,j,nz+1) ; h2 = e(i+1,j,nkmb+1)-e(i+1,j,nz+1)
881  if (e(i,j,nz+1) < e(i+1,j,nz+1)) then
882  if (h1 > h2) h1 = max(e(i,j,nkmb+1)-e(i+1,j,nz+1), h2)
883  elseif (e(i+1,j,nz+1) < e(i,j,nz+1)) then
884  if (h2 > h1) h2 = max(e(i+1,j,nkmb+1)-e(i,j,nz+1), h1)
885  endif
886  h_def_u(i,j) = 0.5*(h1-h2)**2 / ((h1 + h2) + h_neglect)
887  h_norm_u(i,j) = 0.5*(h1+h2)
888  enddo ; enddo
889  if (present(def_rat_u_2lay)) then ; do j=js,je ; do i=is-1,ie
890  ! This is a particular metric of the aggregation into two layers.
891  h1 = e(i,j,1)-e(i,j,nkmb+1) ; h2 = e(i+1,j,1)-e(i+1,j,nkmb+1)
892  if (e(i,j,nkmb+1) < e(i+1,j,nz+1)) then
893  if (h1 > h2) h1 = max(e(i,j,1)-e(i+1,j,nz+1), h2)
894  elseif (e(i+1,j,nkmb+1) < e(i,j,nz+1)) then
895  if (h2 > h1) h2 = max(e(i+1,j,1)-e(i,j,nz+1), h1)
896  endif
897  h_def2_u(i,j) = h_def_u(i,j) + 0.5*(h1-h2)**2 / ((h1 + h2) + h_neglect)
898  enddo ; enddo ; endif
899  do k=1,nkmb ; do j=js,je ; do i=is-1,ie
900  if (present(h)) then
901  h1 = h(i,j,k) ; h2 = h(i+1,j,k)
902  else
903  h1 = e(i,j,k)-e(i,j,k+1) ; h2 = e(i+1,j,k)-e(i+1,j,k+1)
904  endif
905  ! Thickness deficits can not arise simply because a layer's bottom is bounded
906  ! by the bathymetry.
907  if (e(i,j,k+1) < e(i+1,j,nz+1)) then
908  if (h1 > h2) h1 = max(e(i,j,k)-e(i+1,j,nz+1), h2)
909  elseif (e(i+1,j,k+1) < e(i,j,nz+1)) then
910  if (h2 > h1) h2 = max(e(i+1,j,k)-e(i,j,nz+1), h1)
911  endif
912  h_def_u(i,j) = h_def_u(i,j) + 0.5*(h1-h2)**2 / ((h1 + h2) + h_neglect)
913  h_norm_u(i,j) = h_norm_u(i,j) + 0.5*(h1+h2)
914  enddo ; enddo ; enddo
915  if (present(def_rat_u_2lay)) then ; do j=js,je ; do i=is-1,ie
916  def_rat_u(i,j) = g%mask2dCu(i,j) * h_def_u(i,j) / &
917  (max(cs%Hmix_min, h_norm_u(i,j)) + h_neglect)
918  def_rat_u_2lay(i,j) = g%mask2dCu(i,j) * h_def2_u(i,j) / &
919  (max(cs%Hmix_min, h_norm_u(i,j)) + h_neglect)
920  enddo ; enddo ; else ; do j=js,je ; do i=is-1,ie
921  def_rat_u(i,j) = g%mask2dCu(i,j) * h_def_u(i,j) / &
922  (max(cs%Hmix_min, h_norm_u(i,j)) + h_neglect)
923  enddo ; enddo ; endif
924 
925  ! Determine which meridional faces are problematic.
926  do j=js-1,je ; do i=is,ie
927  ! Aggregate all water below the mixed and buffer layers for the purposes of
928  ! this diagnostic.
929  h1 = e(i,j,nkmb+1)-e(i,j,nz+1) ; h2 = e(i,j+1,nkmb+1)-e(i,j+1,nz+1)
930  if (e(i,j,nz+1) < e(i,j+1,nz+1)) then
931  if (h1 > h2) h1 = max(e(i,j,nkmb+1)-e(i,j+1,nz+1), h2)
932  elseif (e(i,j+1,nz+1) < e(i,j,nz+1)) then
933  if (h2 > h1) h2 = max(e(i,j+1,nkmb+1)-e(i,j,nz+1), h1)
934  endif
935  h_def_v(i,j) = 0.5*(h1-h2)**2 / ((h1 + h2) + h_neglect)
936  h_norm_v(i,j) = 0.5*(h1+h2)
937  enddo ; enddo
938  if (present(def_rat_v_2lay)) then ; do j=js-1,je ; do i=is,ie
939  ! This is a particular metric of the aggregation into two layers.
940  h1 = e(i,j,1)-e(i,j,nkmb+1) ; h2 = e(i,j+1,1)-e(i,j+1,nkmb+1)
941  if (e(i,j,nkmb+1) < e(i,j+1,nz+1)) then
942  if (h1 > h2) h1 = max(e(i,j,1)-e(i,j+1,nz+1), h2)
943  elseif (e(i,j+1,nkmb+1) < e(i,j,nz+1)) then
944  if (h2 > h1) h2 = max(e(i,j+1,1)-e(i,j,nz+1), h1)
945  endif
946  h_def2_v(i,j) = h_def_v(i,j) + 0.5*(h1-h2)**2 / ((h1 + h2) + h_neglect)
947  enddo ; enddo ; endif
948  do k=1,nkmb ; do j=js-1,je ; do i=is,ie
949  if (present(h)) then
950  h1 = h(i,j,k) ; h2 = h(i,j+1,k)
951  else
952  h1 = e(i,j,k)-e(i,j,k+1) ; h2 = e(i,j+1,k)-e(i,j+1,k+1)
953  endif
954  ! Thickness deficits can not arise simply because a layer's bottom is bounded
955  ! by the bathymetry.
956  if (e(i,j,k+1) < e(i,j+1,nz+1)) then
957  if (h1 > h2) h1 = max(e(i,j,k)-e(i,j+1,nz+1), h2)
958  elseif (e(i,j+1,k+1) < e(i,j,nz+1)) then
959  if (h2 > h1) h2 = max(e(i,j+1,k)-e(i,j,nz+1), h1)
960  endif
961  h_def_v(i,j) = h_def_v(i,j) + 0.5*(h1-h2)**2 / ((h1 + h2) + h_neglect)
962  h_norm_v(i,j) = h_norm_v(i,j) + 0.5*(h1+h2)
963  enddo ; enddo ; enddo
964  if (present(def_rat_v_2lay)) then ; do j=js-1,je ; do i=is,ie
965  def_rat_v(i,j) = g%mask2dCv(i,j) * h_def_v(i,j) / &
966  (max(cs%Hmix_min, h_norm_v(i,j)) + h_neglect)
967  def_rat_v_2lay(i,j) = g%mask2dCv(i,j) * h_def2_v(i,j) / &
968  (max(cs%Hmix_min, h_norm_v(i,j)) + h_neglect)
969  enddo ; enddo ; else ; do j=js-1,je ; do i=is,ie
970  def_rat_v(i,j) = g%mask2dCv(i,j) * h_def_v(i,j) / &
971  (max(cs%Hmix_min, h_norm_v(i,j)) + h_neglect)
972  enddo ; enddo ; endif
973 
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:19
Here is the caller graph for this function:

◆ 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(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, in H (usually m or kg m-2).
[in,out]tvA structure containing pointers to any available thermodynamic fields. Absent fields have NULL ptrs.
[in]dtTime increment, in s.
[in,out]eaThe amount of fluid moved downward into a
[in,out]ebThe amount of fluid moved upward into a layer;
csThe control structure returned by a previous call to regularize_layers_init.

Definition at line 118 of file MOM_regularize_layers.F90.

References id_clock_pass, mom_error_handler::mom_error(), and regularize_surface().

Referenced by mom_diabatic_driver::diabatic().

118  type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.
119  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
120  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
121  intent(inout) :: h !< Layer thicknesses, in H (usually m or kg m-2).
122  type(thermo_var_ptrs), intent(inout) :: tv !< A structure containing pointers to any
123  !! available thermodynamic fields. Absent fields
124  !! have NULL ptrs.
125  real, intent(in) :: dt !< Time increment, in s.
126  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
127  intent(inout) :: ea !< The amount of fluid moved downward into a
128  !! layer; this should be increased due to mixed
129  !! layer detrainment, in the same units as
130  !! h - usually m or kg m-2 (i.e., H).
131  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
132  intent(inout) :: eb !< The amount of fluid moved upward into a layer;
133  !! this should be increased due to mixed layer
134  !! entrainment, in the same units as h - usually
135  !! m or kg m-2 (i.e., H).
136  type(regularize_layers_cs), pointer :: cs !< The control structure returned by a previous
137  !! call to regularize_layers_init.
138 
139 ! This subroutine partially steps the bulk mixed layer model.
140 ! The following processes are executed, in the order listed.
141 
142 ! Arguments: h - Layer thickness, in m or kg m-2. (Intent in/out)
143 ! The units of h are referred to as H below.
144 ! (in/out) tv - A structure containing pointers to any available
145 ! thermodynamic fields. Absent fields have NULL ptrs.
146 ! (in) dt - Time increment, in s.
147 ! (in/out) ea - The amount of fluid moved downward into a layer; this should
148 ! be increased due to mixed layer detrainment, in the same units
149 ! as h - usually m or kg m-2 (i.e., H).
150 ! (in/out) eb - The amount of fluid moved upward into a layer; this should
151 ! be increased due to mixed layer entrainment, in the same units
152 ! as h - usually m or kg m-2 (i.e., H).
153 ! (in) G - The ocean's grid structure.
154 ! (in) GV - The ocean's vertical grid structure.
155 ! (in) CS - The control structure returned by a previous call to
156 ! regularize_layers_init.
157 
158  integer :: i, j, k, is, ie, js, je, nz
159 
160  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
161 
162  if (.not. associated(cs)) call mom_error(fatal, "MOM_regularize_layers: "//&
163  "Module must be initialized before it is used.")
164 
165  if (cs%regularize_surface_layers) then
166  call cpu_clock_begin(id_clock_pass)
167  call create_group_pass(cs%pass_h,h,g%Domain)
168  call do_group_pass(cs%pass_h,g%Domain)
169  call cpu_clock_end(id_clock_pass)
170  endif
171 
172  if (cs%regularize_surface_layers) then
173  call regularize_surface(h, tv, dt, ea, eb, g, gv, cs)
174  endif
175 
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:19
Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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(param_file_type), intent(in)  param_file,
type(diag_ctrl), intent(inout), target  diag,
type(regularize_layers_cs), pointer  CS 
)
Parameters
[in]timeThe current model time.
[in]gThe ocean's 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 977 of file MOM_regularize_layers.F90.

References id_clock_eos, id_clock_pass, and mom_error_handler::mom_error().

Referenced by mom_diabatic_driver::diabatic_driver_init().

977  type(time_type), target, intent(in) :: time !< The current model time.
978  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
979  type(param_file_type), intent(in) :: param_file !< A structure to parse for
980  !! run-time parameters.
981  type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate
982  !! diagnostic output.
983  type(regularize_layers_cs), pointer :: cs !< A pointer that is set to point to the
984  !! control structure for this module.
985 ! Arguments: Time - The current model time.
986 ! (in) G - The ocean's grid structure.
987 ! (in) param_file - A structure indicating the open file to parse for
988 ! model parameter values.
989 ! (in) diag - A structure that is used to regulate diagnostic output.
990 ! (in/out) CS - A pointer that is set to point to the control structure
991 ! for this module
992 ! This include declares and sets the variable "version".
993 #include "version_variable.h"
994  character(len=40) :: mdl = "MOM_regularize_layers" ! This module's name.
995  logical :: use_temperature
996  integer :: isd, ied, jsd, jed
997  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
998 
999  if (associated(cs)) then
1000  call mom_error(warning, "regularize_layers_init called with an associated"// &
1001  "associated control structure.")
1002  return
1003  else ; allocate(cs) ; endif
1004 
1005  cs%diag => diag
1006  cs%Time => time
1007 
1008 ! Set default, read and log parameters
1009  call log_version(param_file, mdl, version, "")
1010  call get_param(param_file, mdl, "REGULARIZE_SURFACE_LAYERS", cs%regularize_surface_layers, &
1011  "If defined, vertically restructure the near-surface \n"//&
1012  "layers when they have too much lateral variations to \n"//&
1013  "allow for sensible lateral barotropic transports.", &
1014  default=.false.)
1015  if (cs%regularize_surface_layers) then
1016  call get_param(param_file, mdl, "REGULARIZE_SURFACE_DETRAIN", cs%reg_sfc_detrain, &
1017  "If true, allow the buffer layers to detrain into the \n"//&
1018  "interior as a part of the restructuring when \n"//&
1019  "REGULARIZE_SURFACE_LAYERS is true.", default=.true.)
1020  endif
1021 
1022  call get_param(param_file, mdl, "HMIX_MIN", cs%Hmix_min, &
1023  "The minimum mixed layer depth if the mixed layer depth \n"//&
1024  "is determined dynamically.", units="m", default=0.0)
1025  call get_param(param_file, mdl, "REG_SFC_DEFICIT_TOLERANCE", cs%h_def_tol1, &
1026  "The value of the relative thickness deficit at which \n"//&
1027  "to start modifying the layer structure when \n"//&
1028  "REGULARIZE_SURFACE_LAYERS is true.", units="nondim", &
1029  default=0.5)
1030  cs%h_def_tol2 = 0.2 + 0.8*cs%h_def_tol1
1031  cs%h_def_tol3 = 0.3 + 0.7*cs%h_def_tol1
1032  cs%h_def_tol4 = 0.5 + 0.5*cs%h_def_tol1
1033 
1034  call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false.)
1035 ! if (.not. CS%debug) &
1036 ! call get_param(param_file, mdl, "DEBUG_CONSERVATION", CS%debug, &
1037 ! "If true, monitor conservation and extrema.", default=.false.)
1038 
1039  call get_param(param_file, mdl, "ALLOW_CLOCKS_IN_OMP_LOOPS", &
1040  cs%allow_clocks_in_omp_loops, &
1041  "If true, clocks can be called from inside loops that can \n"//&
1042  "be threaded. To run with multiple threads, set to False.", &
1043  default=.true.)
1044 
1045  cs%id_def_rat = register_diag_field('ocean_model', 'deficit_ratio', diag%axesT1, &
1046  time, 'Max face thickness deficit ratio', 'Nondim')
1047 
1048 #ifdef DEBUG_CODE
1049  cs%id_def_rat_2 = register_diag_field('ocean_model', 'deficit_rat2', diag%axesT1, &
1050  time, 'Corrected thickness deficit ratio', 'Nondim')
1051  cs%id_def_rat_3 = register_diag_field('ocean_model', 'deficit_rat3', diag%axesT1, &
1052  time, 'Filtered thickness deficit ratio', 'Nondim')
1053  cs%id_e1 = register_diag_field('ocean_model', 'er_1', diag%axesTi, &
1054  time, 'Intial interface depths before remapping', 'm')
1055  cs%id_e2 = register_diag_field('ocean_model', 'er_2', diag%axesTi, &
1056  time, 'Intial interface depths after remapping', 'm')
1057  cs%id_e3 = register_diag_field('ocean_model', 'er_3', diag%axesTi, &
1058  time, 'Intial interface depths filtered', 'm')
1059 
1060  cs%id_def_rat_u = register_diag_field('ocean_model', 'defrat_u', diag%axesCu1, &
1061  time, 'U-point thickness deficit ratio', 'Nondim')
1062  cs%id_def_rat_u_1b = register_diag_field('ocean_model', 'defrat_u_1b', diag%axesCu1, &
1063  time, 'U-point 2-layer thickness deficit ratio', 'Nondim')
1064  cs%id_def_rat_u_2 = register_diag_field('ocean_model', 'defrat_u_2', diag%axesCu1, &
1065  time, 'U-point corrected thickness deficit ratio', 'Nondim')
1066  cs%id_def_rat_u_2b = register_diag_field('ocean_model', 'defrat_u_2b', diag%axesCu1, &
1067  time, 'U-point corrected 2-layer thickness deficit ratio', 'Nondim')
1068  cs%id_def_rat_u_3 = register_diag_field('ocean_model', 'defrat_u_3', diag%axesCu1, &
1069  time, 'U-point filtered thickness deficit ratio', 'Nondim')
1070  cs%id_def_rat_u_3b = register_diag_field('ocean_model', 'defrat_u_3b', diag%axesCu1, &
1071  time, 'U-point filtered 2-layer thickness deficit ratio', 'Nondim')
1072 
1073  cs%id_def_rat_v = register_diag_field('ocean_model', 'defrat_v', diag%axesCv1, &
1074  time, 'V-point thickness deficit ratio', 'Nondim')
1075  cs%id_def_rat_v_1b = register_diag_field('ocean_model', 'defrat_v_1b', diag%axesCv1, &
1076  time, 'V-point 2-layer thickness deficit ratio', 'Nondim')
1077  cs%id_def_rat_v_2 = register_diag_field('ocean_model', 'defrat_v_2', diag%axesCv1, &
1078  time, 'V-point corrected thickness deficit ratio', 'Nondim')
1079  cs%id_def_rat_v_2b = register_diag_field('ocean_model', 'defrat_v_2b', diag%axesCv1, &
1080  time, 'V-point corrected 2-layer thickness deficit ratio', 'Nondim')
1081  cs%id_def_rat_v_3 = register_diag_field('ocean_model', 'defrat_v_3', diag%axesCv1, &
1082  time, 'V-point filtered thickness deficit ratio', 'Nondim')
1083  cs%id_def_rat_v_3b = register_diag_field('ocean_model', 'defrat_v_3b', diag%axesCv1, &
1084  time, 'V-point filtered 2-layer thickness deficit ratio', 'Nondim')
1085 #endif
1086 
1087  if(cs%allow_clocks_in_omp_loops) then
1088  id_clock_eos = cpu_clock_id('(Ocean regularize_layers EOS)', grain=clock_routine)
1089  endif
1090  id_clock_pass = cpu_clock_id('(Ocean regularize_layers halo updates)', grain=clock_routine)
1091 
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:19
Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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(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, in H (usually m or kg m-2).
[in,out]tvA structure containing pointers to any available thermodynamic fields. Absent fields have NULL ptrs.
[in]dtTime increment, in s.
[in,out]eaThe amount of fluid moved downward into a
[in,out]ebThe amount of fluid moved upward into a layer;
csThe control structure returned by a previous call to regularize_layers_init.

Definition at line 181 of file MOM_regularize_layers.F90.

References find_deficit_ratios(), id_clock_eos, and mom_error_handler::mom_error().

Referenced by regularize_layers().

181  type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.
182  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
183  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
184  intent(inout) :: h !< Layer thicknesses, in H (usually m or kg m-2).
185  type(thermo_var_ptrs), intent(inout) :: tv !< A structure containing pointers to any
186  !! available thermodynamic fields. Absent fields
187  !! have NULL ptrs.
188  real, intent(in) :: dt !< Time increment, in s.
189  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
190  intent(inout) :: ea !< The amount of fluid moved downward into a
191  !! layer; this should be increased due to mixed
192  !! layer detrainment, in the same units as h -
193  !! usually m or kg m-2 (i.e., H).
194  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
195  intent(inout) :: eb !< The amount of fluid moved upward into a layer;
196  !! this should be increased due to mixed layer
197  !! entrainment, in the same units as h - usually
198  !! m or kg m-2 (i.e., H).
199  type(regularize_layers_cs), pointer :: cs !< The control structure returned by a previous
200  !! call to regularize_layers_init.
201 
202 ! This subroutine ensures that there is a degree of horizontal smoothness
203 ! in the depths of the near-surface interfaces.
204 
205 ! Arguments: h - Layer thickness, in m or kg m-2. (Intent in/out)
206 ! The units of h are referred to as H below.
207 ! (in/out) tv - A structure containing pointers to any available
208 ! thermodynamic fields. Absent fields have NULL ptrs.
209 ! (in) dt - Time increment, in s.
210 ! (in/out) ea - The amount of fluid moved downward into a layer; this should
211 ! be increased due to mixed layer detrainment, in the same units
212 ! as h - usually m or kg m-2 (i.e., H).
213 ! (in/out) eb - The amount of fluid moved upward into a layer; this should
214 ! be increased due to mixed layer entrainment, in the same units
215 ! as h - usually m or kg m-2 (i.e., H).
216 ! (in) G - The ocean's grid structure.
217 ! (in) GV - The ocean's vertical grid structure.
218 ! (in) CS - The control structure returned by a previous call to
219 ! regularize_layers_init.
220 
221  real, dimension(SZIB_(G),SZJ_(G)) :: &
222  def_rat_u ! The ratio of the thickness deficit to the minimum depth, ND.
223  real, dimension(SZI_(G),SZJB_(G)) :: &
224  def_rat_v ! The ratio of the thickness deficit to the minimum depth, ND.
225  real, dimension(SZI_(G),SZJ_(G)) :: &
226  def_rat_h ! The ratio of the thickness deficit to the minimum depth, ND.
227  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: &
228  e ! The interface depths, in H, positive upward.
229 
230 #ifdef DEBUG_CODE
231  real, dimension(SZIB_(G),SZJ_(G)) :: &
232  def_rat_u_1b, def_rat_u_2, def_rat_u_2b, def_rat_u_3, def_rat_u_3b
233  real, dimension(SZI_(G),SZJB_(G)) :: &
234  def_rat_v_1b, def_rat_v_2, def_rat_v_2b, def_rat_v_3, def_rat_v_3b
235  real, dimension(SZI_(G),SZJB_(G)) :: &
236  def_rat_h2, def_rat_h3
237  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: &
238  ef ! The filtered interface depths, in H, positive upward.
239 #endif
240 
241  real, dimension(SZI_(G),SZK_(G)+1) :: &
242  e_filt, e_2d ! The interface depths, in H, positive upward.
243  real, dimension(SZI_(G),SZK_(G)) :: &
244  h_2d, & ! A 2-d version of h, in H.
245  t_2d, & ! A 2-d version of tv%T, in deg C.
246  s_2d, & ! A 2-d version of tv%S, in PSU.
247  rcv, & ! A 2-d version of the coordinate density, in kg m-3.
248  h_2d_init, & ! The initial value of h_2d, in H.
249  t_2d_init, & ! THe initial value of T_2d, in deg C.
250  s_2d_init, & ! The initial value of S_2d, in PSU.
251  d_eb, & ! The downward increase across a layer in the entrainment from
252  ! below, in H. The sign convention is that positive values of
253  ! d_eb correspond to a gain in mass by a layer by upward motion.
254  d_ea ! The upward increase across a layer in the entrainment from
255  ! above, in H. The sign convention is that positive values of
256  ! d_ea mean a net gain in mass by a layer from downward motion.
257  real, dimension(SZI_(G)) :: &
258  p_ref_cv, & ! Reference pressure for the potential density which defines
259  ! the coordinate variable, set to P_Ref, in Pa.
260  rcv_tol, & ! A tolerence, relative to the target density differences
261  ! between layers, for detraining into the interior, ND.
262  h_add_tgt, h_add_tot, &
263  h_tot1, th_tot1, sh_tot1, &
264  h_tot3, th_tot3, sh_tot3, &
265  h_tot2, th_tot2, sh_tot2
266  real, dimension(SZK_(G)) :: &
267  h_prev_1d ! The previous thicknesses, in H.
268  real :: i_dtol ! The inverse of the tolerance changes, nondim.
269  real :: i_dtol34 ! The inverse of the tolerance changes, nondim.
270  real :: h1, h2 ! Temporary thicknesses, in H.
271  real :: e_e, e_w, e_n, e_s ! Temporary interface heights, in H.
272  real :: wt ! The weight of the filted interfaces in setting the targets, ND.
273  real :: scale ! A scaling factor, ND.
274  real :: h_neglect ! A thickness that is so small it is usually lost
275  ! in roundoff and can be neglected, in H.
276  real, dimension(SZK_(G)+1) :: &
277  int_flux, int_tflux, int_sflux, int_rflux
278  real :: h_add
279  real :: h_det_tot
280  real :: max_def_rat
281  real :: rcv_min_det ! The lightest (min) and densest (max) coordinate density
282  real :: rcv_max_det ! that can detrain into a layer, in kg m-3.
283 
284  real :: int_top, int_bot
285  real :: h_predicted
286  real :: h_prev
287  real :: h_deficit
288 
289  logical :: cols_left, ent_any, more_ent_i(szi_(g)), ent_i(szi_(g))
290  logical :: det_any, det_i(szi_(g))
291  logical :: do_j(szj_(g)), do_i(szi_(g)), find_i(szi_(g))
292  logical :: debug = .false.
293  logical :: fatal_error
294  character(len=256) :: mesg ! Message for error messages.
295  integer :: i, j, k, is, ie, js, je, nz, nkmb, nkml, k1, k2, k3, ks, nz_filt
296 
297  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
298 
299  if (.not. associated(cs)) call mom_error(fatal, "MOM_regularize_layers: "//&
300  "Module must be initialized before it is used.")
301 
302  if (gv%nkml<1) return
303  nkmb = gv%nk_rho_varies ; nkml = gv%nkml
304  if (.not.ASSOCIATED(tv%eqn_of_state)) call mom_error(fatal, &
305  "MOM_regularize_layers: This module now requires the use of temperature and "//&
306  "an equation of state.")
307 
308  h_neglect = gv%H_subroundoff
309  debug = (debug .or. cs%debug)
310 #ifdef DEBUG_CODE
311  debug = .true.
312  if (cs%id_def_rat_2 > 0) then ! Calculate over a slightly larger domain.
313  is = g%isc-1 ; ie = g%iec+1 ; js = g%jsc-1 ; je = g%jec+1
314  endif
315 #endif
316 
317  i_dtol = 1.0 / max(cs%h_def_tol2 - cs%h_def_tol1, 1e-40)
318  i_dtol34 = 1.0 / max(cs%h_def_tol4 - cs%h_def_tol3, 1e-40)
319 
320  p_ref_cv(:) = tv%P_Ref
321 
322  do j=js-1,je+1 ; do i=is-1,ie+1
323  e(i,j,1) = 0.0
324  enddo ; enddo
325  do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
326  e(i,j,k+1) = e(i,j,k) - h(i,j,k)
327  enddo ; enddo ; enddo
328 
329 #ifdef DEBUG_CODE
330  call find_deficit_ratios(e, def_rat_u, def_rat_v, g, gv, cs, def_rat_u_1b, def_rat_v_1b, 1, h)
331 #else
332  call find_deficit_ratios(e, def_rat_u, def_rat_v, g, gv, cs, h=h)
333 #endif
334  ! Determine which columns are problematic
335  do j=js,je ; do_j(j) = .false. ; enddo
336  do j=js,je ; do i=is,ie
337  def_rat_h(i,j) = max(def_rat_u(i-1,j), def_rat_u(i,j), &
338  def_rat_v(i,j-1), def_rat_v(i,j))
339  if (def_rat_h(i,j) > cs%h_def_tol1) do_j(j) = .true.
340  enddo ; enddo
341 
342 #ifdef DEBUG_CODE
343  if ((cs%id_def_rat_3 > 0) .or. (cs%id_e3 > 0) .or. &
344  (cs%id_def_rat_u_3 > 0) .or. (cs%id_def_rat_u_3b > 0) .or. &
345  (cs%id_def_rat_v_3 > 0) .or. (cs%id_def_rat_v_3b > 0) ) then
346  do j=js-1,je+1 ; do i=is-1,ie+1
347  ef(i,j,1) = 0.0
348  enddo ; enddo
349  do k=2,nz+1 ; do j=js,je ; do i=is,ie
350  if (g%mask2dCu(i,j) <= 0.0) then ; e_e = e(i,j,k) ; else
351  e_e = max(e(i+1,j,k) + min(e(i,j,k) - e(i+1,j,nz+1), 0.0), e(i,j,nz+1))
352  endif
353  if (g%mask2dCu(i-1,j) <= 0.0) then ; e_w = e(i,j,k) ; else
354  e_w = max(e(i-1,j,k) + min(e(i,j,k) - e(i-1,j,nz+1), 0.0), e(i,j,nz+1))
355  endif
356  if (g%mask2dCv(i,j) <= 0.0) then ; e_n = e(i,j,k) ; else
357  e_n = max(e(i,j+1,k) + min(e(i,j,k) - e(i,j+1,nz+1), 0.0), e(i,j,nz+1))
358  endif
359  if (g%mask2dCv(i,j-1) <= 0.0) then ; e_s = e(i,j,k) ; else
360  e_s = max(e(i,j-1,k) + min(e(i,j,k) - e(i,j-1,nz+1), 0.0), e(i,j,nz+1))
361  endif
362 
363  wt = 1.0
364  ef(i,j,k) = (1.0 - 0.5*wt) * e(i,j,k) + &
365  wt * 0.125 * ((e_e + e_w) + (e_n + e_s))
366  enddo ; enddo ; enddo
367  call find_deficit_ratios(ef, def_rat_u_3, def_rat_v_3, g, gv, cs, def_rat_u_3b, def_rat_v_3b)
368 
369  ! Determine which columns are problematic
370  do j=js,je ; do i=is,ie
371  def_rat_h3(i,j) = max(def_rat_u_3(i-1,j), def_rat_u_3(i,j), &
372  def_rat_v_3(i,j-1), def_rat_v_3(i,j))
373  enddo ; enddo
374 
375  if (cs%id_e3 > 0) call post_data(cs%id_e3, ef, cs%diag)
376  if (cs%id_def_rat_3 > 0) call post_data(cs%id_def_rat_3, def_rat_h3, cs%diag)
377  if (cs%id_def_rat_u_3 > 0) call post_data(cs%id_def_rat_u_3, def_rat_u_3, cs%diag)
378  if (cs%id_def_rat_u_3b > 0) call post_data(cs%id_def_rat_u_3b, def_rat_u_3b, cs%diag)
379  if (cs%id_def_rat_v_3 > 0) call post_data(cs%id_def_rat_v_3, def_rat_v_3, cs%diag)
380  if (cs%id_def_rat_v_3b > 0) call post_data(cs%id_def_rat_v_3b, def_rat_v_3b, cs%diag)
381  endif
382 #endif
383 
384 
385  ! Now restructure the layers.
386 !$OMP parallel do default(none) shared(is,ie,js,je,nz,do_j,def_rat_h,CS,nkmb,G,GV,&
387 !$OMP e,I_dtol,h,tv,debug,h_neglect,p_ref_cv,ea, &
388 !$OMP eb,id_clock_EOS,nkml) &
389 !$OMP private(d_ea,d_eb,max_def_rat,do_i,nz_filt,e_e,e_w,&
390 !$OMP e_n,e_s,wt,e_filt,e_2d,h_2d,T_2d,S_2d, &
391 !$OMP h_2d_init,T_2d_init,S_2d_init,ent_any, &
392 !$OMP more_ent_i,ent_i,h_add_tgt,h_add_tot, &
393 !$OMP cols_left,h_add,h_prev,ks,det_any,det_i, &
394 !$OMP Rcv_tol,Rcv,k1,k2,h_det_tot,Rcv_min_det, &
395 !$OMP Rcv_max_det,h_deficit,h_tot3,Th_tot3, &
396 !$OMP Sh_tot3,scale,int_top,int_flux,int_Rflux, &
397 !$OMP int_Tflux,int_Sflux,int_bot,h_prev_1d, &
398 !$OMP h_tot1,Th_tot1,Sh_tot1,h_tot2,Th_tot2, &
399 !$OMP Sh_tot2,h_predicted,fatal_error,mesg )
400  do j=js,je ; if (do_j(j)) then
401 
402 ! call cpu_clock_begin(id_clock_EOS)
403 ! call calculate_density_derivs(T(:,1), S(:,1), p_ref_cv, dRcv_dT, dRcv_dS, &
404 ! is, ie-is+1, tv%eqn_of_state)
405 ! call cpu_clock_end(id_clock_EOS)
406 
407  do k=1,nz ; do i=is,ie ; d_ea(i,k) = 0.0 ; d_eb(i,k) = 0.0 ; enddo ; enddo
408 
409  max_def_rat = 0.0
410  do i=is,ie
411  do_i(i) = def_rat_h(i,j) > cs%h_def_tol1
412  if (def_rat_h(i,j) > max_def_rat) max_def_rat = def_rat_h(i,j)
413  enddo
414  nz_filt = nkmb+1 ; if (max_def_rat > cs%h_def_tol3) nz_filt = nz+1
415 
416  ! Find a 2-D 1-2-1 filtered version of e to target. Area weights are
417  ! deliberately omitted here. This is slightly more complicated than a
418  ! simple filter so that the effects of topography are eliminated.
419  do k=1,nz_filt ; do i=is,ie ; if (do_i(i)) then
420  if (g%mask2dCu(i,j) <= 0.0) then ; e_e = e(i,j,k) ; else
421  e_e = max(e(i+1,j,k) + min(e(i,j,k) - e(i+1,j,nz+1), 0.0), &
422  e(i,j,nz+1) + (nz+1-k)*gv%Angstrom)
423 
424  endif
425  if (g%mask2dCu(i-1,j) <= 0.0) then ; e_w = e(i,j,k) ; else
426  e_w = max(e(i-1,j,k) + min(e(i,j,k) - e(i-1,j,nz+1), 0.0), &
427  e(i,j,nz+1) + (nz+1-k)*gv%Angstrom)
428  endif
429  if (g%mask2dCv(i,j) <= 0.0) then ; e_n = e(i,j,k) ; else
430  e_n = max(e(i,j+1,k) + min(e(i,j,k) - e(i,j+1,nz+1), 0.0), &
431  e(i,j,nz+1) + (nz+1-k)*gv%Angstrom)
432  endif
433  if (g%mask2dCv(i,j-1) <= 0.0) then ; e_s = e(i,j,k) ; else
434  e_s = max(e(i,j-1,k) + min(e(i,j,k) - e(i,j-1,nz+1), 0.0), &
435  e(i,j,nz+1) + (nz+1-k)*gv%Angstrom)
436  endif
437 
438  wt = max(0.0, min(1.0, i_dtol*(def_rat_h(i,j)-cs%h_def_tol1)))
439 
440  e_filt(i,k) = (1.0 - 0.5*wt) * e(i,j,k) + &
441  wt * 0.125 * ((e_e + e_w) + (e_n + e_s))
442  e_2d(i,k) = e(i,j,k)
443  endif ; enddo ; enddo
444  do k=1,nz ; do i=is,ie
445  h_2d(i,k) = h(i,j,k)
446  t_2d(i,k) = tv%T(i,j,k) ; s_2d(i,k) = tv%S(i,j,k)
447  enddo ; enddo
448 
449  if (debug) then
450  do k=1,nz ; do i=is,ie ; if (do_i(i)) then
451  h_2d_init(i,k) = h(i,j,k)
452  t_2d_init(i,k) = tv%T(i,j,k) ; s_2d_init(i,k) = tv%S(i,j,k)
453  endif ; enddo ; enddo
454  endif
455 
456  ! First, try to entrain from the interior.
457  ent_any = .false.
458  do i=is,ie
459  more_ent_i(i) = .false. ; ent_i(i) = .false.
460  h_add_tgt(i) = 0.0 ; h_add_tot(i) = 0.0
461  if (do_i(i) .and. (e_2d(i,nkmb+1) > e_filt(i,nkmb+1))) then
462  more_ent_i(i) = .true. ; ent_i(i) = .true. ; ent_any = .true.
463  h_add_tgt(i) = e_2d(i,nkmb+1) - e_filt(i,nkmb+1)
464  endif
465  enddo
466 
467  if (ent_any) then
468  do k=nkmb+1,nz
469  cols_left = .false.
470  do i=is,ie ; if (more_ent_i(i)) then
471  if (h_2d(i,k) - gv%Angstrom > h_neglect) then
472  if (e_2d(i,nkmb+1)-e_filt(i,nkmb+1) > h_2d(i,k) - gv%Angstrom) then
473  h_add = h_2d(i,k) - gv%Angstrom
474  h_2d(i,k) = gv%Angstrom
475  else
476  h_add = e_2d(i,nkmb+1)-e_filt(i,nkmb+1)
477  h_2d(i,k) = h_2d(i,k) - h_add
478  endif
479  d_eb(i,k-1) = d_eb(i,k-1) + h_add
480  h_add_tot(i) = h_add_tot(i) + h_add
481  e_2d(i,nkmb+1) = e_2d(i,nkmb+1) - h_add
482  h_prev = h_2d(i,nkmb)
483  h_2d(i,nkmb) = h_2d(i,nkmb) + h_add
484 
485  t_2d(i,nkmb) = (h_prev*t_2d(i,nkmb) + h_add*t_2d(i,k)) / h_2d(i,nkmb)
486  s_2d(i,nkmb) = (h_prev*s_2d(i,nkmb) + h_add*s_2d(i,k)) / h_2d(i,nkmb)
487 
488  if ((e_2d(i,nkmb+1) <= e_filt(i,nkmb+1)) .or. &
489  (h_add_tot(i) > 0.6*h_add_tgt(i))) then !### 0.6 is adjustable?.
490  more_ent_i(i) = .false.
491  else
492  cols_left = .true.
493  endif
494  else
495  cols_left = .true.
496  endif
497  endif ; enddo
498  if (.not.cols_left) exit
499  enddo
500 
501  ks = min(k-1,nz-1)
502  do k=ks,nkmb,-1 ; do i=is,ie ; if (ent_i(i)) then
503  d_eb(i,k) = d_eb(i,k) + d_eb(i,k+1)
504  endif ; enddo ; enddo
505  endif ! ent_any
506 
507  ! This is where code to detrain to the interior will go.
508  ! The buffer layers can only detrain water into layers when the buffer
509  ! layer potential density is between (c*Rlay(k-1) + (1-c)*Rlay(k)) and
510  ! (c*Rlay(k+1) + (1-c)*Rlay(k)), where 0.5 <= c < 1.0.
511  ! Do not detrain if the 2-layer deficit ratio is not significant.
512  ! Detrainment must be able to come from all mixed and buffer layers.
513  ! All water is moved out of the buffer layers below before moving from
514  ! a shallower layer (characteristics do not cross).
515  det_any = .false.
516  if ((max_def_rat > cs%h_def_tol3) .and. (cs%reg_sfc_detrain)) then
517  do i=is,ie
518  det_i(i) = .false. ; rcv_tol(i) = 0.0
519  if (do_i(i) .and. (e_2d(i,nkmb+1) < e_filt(i,nkmb+1)) .and. &
520  (def_rat_h(i,j) > cs%h_def_tol3)) then
521  det_i(i) = .true. ; det_any = .true.
522  rcv_tol(i) = min((def_rat_h(i,j) - cs%h_def_tol3), 1.0)
523  endif
524  enddo
525  endif
526  if (det_any) then
527  call cpu_clock_begin(id_clock_eos)
528  do k=1,nkmb
529  call calculate_density(t_2d(:,k),s_2d(:,k),p_ref_cv,rcv(:,k), &
530  is,ie-is+1,tv%eqn_of_state)
531  enddo
532  call cpu_clock_end(id_clock_eos)
533 
534  do i=is,ie ; if (det_i(i)) then
535  k1 = nkmb ; k2 = nz
536  h_det_tot = 0.0
537  do ! This loop is terminated by exits.
538  if (k1 <= 1) exit
539  if (k2 <= nkmb) exit
540  ! ### The 0.6 here should be adjustable? It gives 20% overlap for now.
541  rcv_min_det = gv%Rlay(k2) + 0.6*rcv_tol(i)*(gv%Rlay(k2-1)-gv%Rlay(k2))
542  if (k2 < nz) then
543  rcv_max_det = gv%Rlay(k2) + 0.6*rcv_tol(i)*(gv%Rlay(k2+1)-gv%Rlay(k2))
544  else
545  rcv_max_det = gv%Rlay(nz) + 0.6*rcv_tol(i)*(gv%Rlay(nz)-gv%Rlay(nz-1))
546  endif
547  if (rcv(i,k1) > rcv_max_det) &
548  exit ! All shallower interior layers are too light for detrainment.
549 
550  h_deficit = (e_filt(i,k2)-e_filt(i,k2+1)) - h_2d(i,k2)
551  if ((e_filt(i,k2) > e_2d(i,k1+1)) .and. (h_deficit > 0.0) .and. &
552  (rcv(i,k1) < rcv_max_det) .and. (rcv(i,k1) > rcv_min_det)) then
553  ! Detrainment will occur.
554  h_add = min(e_filt(i,k2) - e_2d(i,k2), h_deficit )
555  if (h_add < h_2d(i,k1)) then
556  ! Only part of layer k1 detrains.
557  if (h_add > 0.0) then
558  h_prev = h_2d(i,k2)
559  h_2d(i,k2) = h_2d(i,k2) + h_add
560  e_2d(i,k2) = e_2d(i,k2+1) + h_2d(i,k2)
561  d_ea(i,k2) = d_ea(i,k2) + h_add
562  ! ### THIS IS UPWIND. IT SHOULD BE HIGHER ORDER...
563  t_2d(i,k2) = (h_prev*t_2d(i,k2) + h_add*t_2d(i,k1)) / h_2d(i,k2)
564  s_2d(i,k2) = (h_prev*s_2d(i,k2) + h_add*s_2d(i,k1)) / h_2d(i,k2)
565  h_det_tot = h_det_tot + h_add
566 
567  h_2d(i,k1) = h_2d(i,k1) - h_add
568  do k3=k1,nkmb ; e_2d(i,k3+1) = e_2d(i,k3) - h_2d(i,k3) ; enddo
569  do k3=k1+1,nkmb ; d_ea(i,k3) = d_ea(i,k3) + h_add ; enddo
570  else
571  if (h_add < 0.0) &
572  call mom_error(fatal, "h_add is negative. Some logic is wrong.")
573  h_add = 0.0 ! This usually should not happen...
574  endif
575 
576  ! Move up to the next target layer.
577  k2 = k2-1
578  if (k2>nkmb+1) e_2d(i,k2) = e_2d(i,k2) + h_det_tot
579  else
580  h_add = h_2d(i,k1)
581  h_prev = h_2d(i,k2)
582  h_2d(i,k2) = h_2d(i,k2) + h_add
583  e_2d(i,k2) = e_2d(i,k2+1) + h_2d(i,k2)
584  d_ea(i,k2) = d_ea(i,k2) + h_add
585  t_2d(i,k2) = (h_prev*t_2d(i,k2) + h_add*t_2d(i,k1)) / h_2d(i,k2)
586  s_2d(i,k2) = (h_prev*s_2d(i,k2) + h_add*s_2d(i,k1)) / h_2d(i,k2)
587  h_det_tot = h_det_tot + h_add
588 
589  h_2d(i,k1) = 0.0
590  do k3=k1,nkmb ; e_2d(i,k3+1) = e_2d(i,k3) - h_2d(i,k3) ; enddo
591  do k3=k1+1,nkmb ; d_ea(i,k3) = d_ea(i,k3) + h_add ; enddo
592 
593  ! Move up to the next source layer.
594  k1 = k1-1
595  endif
596 
597  else
598  ! Move up to the next target layer.
599  k2 = k2-1
600  if (k2>nkmb+1) e_2d(i,k2) = e_2d(i,k2) + h_det_tot
601  endif
602 
603  enddo ! exit terminated loop.
604  endif ; enddo
605  ! ### This could be faster if the deepest k with nonzero d_ea were kept.
606  do k=nz-1,nkmb+1,-1 ; do i=is,ie ; if (det_i(i)) then
607  d_ea(i,k) = d_ea(i,k) + d_ea(i,k+1)
608  endif ; enddo ; enddo
609  endif ! Detrainment to the interior.
610  if (debug) then
611  do i=is,ie ; h_tot3(i) = 0.0 ; th_tot3(i) = 0.0 ; sh_tot3(i) = 0.0 ; enddo
612  do k=1,nz ; do i=is,ie ; if (do_i(i)) then
613  h_tot3(i) = h_tot3(i) + h_2d(i,k)
614  th_tot3(i) = th_tot3(i) + h_2d(i,k) * t_2d(i,k)
615  sh_tot3(i) = sh_tot3(i) + h_2d(i,k) * s_2d(i,k)
616  endif ; enddo ; enddo
617  endif
618 
619  do i=is,ie ; if (do_i(i)) then
620  ! Rescale the interface targets so the depth at the bottom of the deepest
621  ! buffer layer matches.
622  scale = e_2d(i,nkmb+1) / e_filt(i,nkmb+1)
623  do k=2,nkmb+1 ; e_filt(i,k) = e_filt(i,k) * scale ; enddo
624 
625  ! Ensure that layer 1 only has water from layers 1 to nkml and rescale
626  ! the remaining layer thicknesses if necessary.
627  if (e_filt(i,2) < e_2d(i,nkml)) then
628  scale = (e_2d(i,nkml) - e_filt(i,nkmb+1)) / &
629  ((e_filt(i,2) - e_filt(i,nkmb+1)) + h_neglect)
630  do k=3,nkmb
631  e_filt(i,k) = e_filt(i,nkmb+1) + scale * (e_filt(i,k) - e_filt(i,nkmb+1))
632  enddo
633  e_filt(i,2) = e_2d(i,nkml)
634  endif
635 
636  ! Map the water back into the layers.
637  k1 = 1 ; k2 = 1
638  int_top = 0.0
639  do k=1,nkmb+1
640  int_flux(k) = 0.0 ; int_rflux(k) = 0.0
641  int_tflux(k) = 0.0 ; int_sflux(k) = 0.0
642  enddo
643  do k=1,2*nkmb
644  int_bot = max(e_2d(i,k1+1),e_filt(i,k2+1))
645  h_add = int_top - int_bot
646 
647  if (k2 > k1) then
648  do k3=k1+1,k2
649  d_ea(i,k3) = d_ea(i,k3) + h_add
650  int_flux(k3) = int_flux(k3) + h_add
651  int_tflux(k3) = int_tflux(k3) + h_add*t_2d(i,k1)
652  int_sflux(k3) = int_sflux(k3) + h_add*s_2d(i,k1)
653  enddo
654  elseif (k1 > k2) then
655  do k3=k2,k1-1
656  d_eb(i,k3) = d_eb(i,k3) + h_add
657  int_flux(k3+1) = int_flux(k3+1) - h_add
658  int_tflux(k3+1) = int_tflux(k3+1) - h_add*t_2d(i,k1)
659  int_sflux(k3+1) = int_sflux(k3+1) - h_add*s_2d(i,k1)
660  enddo
661  endif
662 
663  if (int_bot <= e_filt(i,k2+1)) then
664  ! Increment the target layer.
665  k2 = k2 + 1
666  elseif (int_bot <= e_2d(i,k1+1)) then
667  ! Increment the source layer.
668  k1 = k1 + 1
669  else
670  call mom_error(fatal, &
671  "Regularize_surface: Could not increment target or source.")
672  endif
673  if ((k1 > nkmb) .or. (k2 > nkmb)) exit
674  int_top = int_bot
675  enddo
676  if (k2 < nkmb) &
677  call mom_error(fatal, "Regularize_surface: Did not assign fluid to layer nkmb.")
678 
679  ! Note that movement of water across the base of the bottommost buffer
680  ! layer has already been dealt with separately.
681  do k=1,nkmb ; h_prev_1d(k) = h_2d(i,k) ; enddo
682  h_2d(i,1) = h_2d(i,1) - int_flux(2)
683  do k=2,nkmb-1
684  h_2d(i,k) = h_2d(i,k) + (int_flux(k) - int_flux(k+1))
685  enddo
686  ! Note that movement of water across the base of the bottommost buffer
687  ! layer has already been dealt with separately.
688  h_2d(i,nkmb) = h_2d(i,nkmb) + int_flux(nkmb)
689 
690  t_2d(i,1) = (t_2d(i,1)*h_prev_1d(1) - int_tflux(2)) / h_2d(i,1)
691  s_2d(i,1) = (s_2d(i,1)*h_prev_1d(1) - int_sflux(2)) / h_2d(i,1)
692  do k=2,nkmb-1
693  t_2d(i,k) = (t_2d(i,k)*h_prev_1d(k) + (int_tflux(k) - int_tflux(k+1))) / h_2d(i,k)
694  s_2d(i,k) = (s_2d(i,k)*h_prev_1d(k) + (int_sflux(k) - int_sflux(k+1))) / h_2d(i,k)
695  enddo
696  t_2d(i,nkmb) = (t_2d(i,nkmb)*h_prev_1d(nkmb) + int_tflux(nkmb) ) / h_2d(i,nkmb)
697  s_2d(i,nkmb) = (s_2d(i,nkmb)*h_prev_1d(nkmb) + int_sflux(nkmb) ) / h_2d(i,nkmb)
698 
699  endif ; enddo ! i-loop
700 
701  ! Copy the interior thicknesses and other fields back to the 3-d arrays.
702  do k=1,nz ; do i=is,ie ; if (do_i(i)) then
703  h(i,j,k) = h_2d(i,k)
704  tv%T(i,j,k) = t_2d(i,k) ; tv%S(i,j,k) = s_2d(i,k)
705  ea(i,j,k) = ea(i,j,k) + d_ea(i,k)
706  eb(i,j,k) = eb(i,j,k) + d_eb(i,k)
707  endif ; enddo ; enddo
708 
709  if (debug) then
710  do i=is,ie ; h_tot1(i) = 0.0 ; th_tot1(i) = 0.0 ; sh_tot1(i) = 0.0 ; enddo
711  do i=is,ie ; h_tot2(i) = 0.0 ; th_tot2(i) = 0.0 ; sh_tot2(i) = 0.0 ; enddo
712 
713  do k=1,nz ; do i=is,ie ; if (do_i(i)) then
714  h_tot1(i) = h_tot1(i) + h_2d_init(i,k)
715  h_tot2(i) = h_tot2(i) + h(i,j,k)
716 
717  th_tot1(i) = th_tot1(i) + h_2d_init(i,k) * t_2d_init(i,k)
718  th_tot2(i) = th_tot2(i) + h(i,j,k) * tv%T(i,j,k)
719  sh_tot1(i) = sh_tot1(i) + h_2d_init(i,k) * s_2d_init(i,k)
720  sh_tot2(i) = sh_tot2(i) + h(i,j,k) * tv%S(i,j,k)
721  if (h(i,j,k) < 0.0) &
722  call mom_error(fatal,"regularize_surface: Negative thicknesses.")
723  if (k==1) then ; h_predicted = h_2d_init(i,k) + (d_eb(i,k) - d_ea(i,k+1))
724  elseif (k==nz) then ; h_predicted = h_2d_init(i,k) + (d_ea(i,k) - d_eb(i,k-1))
725  else
726  h_predicted = h_2d_init(i,k) + ((d_ea(i,k) - d_eb(i,k-1)) + &
727  (d_eb(i,k) - d_ea(i,k+1)))
728  endif
729  if (abs(h(i,j,k) - h_predicted) > max(1e-9*abs(h_predicted),gv%Angstrom)) &
730  call mom_error(fatal, "regularize_surface: d_ea mismatch.")
731  endif ; enddo ; enddo
732  do i=is,ie ; if (do_i(i)) then
733  fatal_error = .false.
734  if (abs(h_tot1(i) - h_tot2(i)) > 1e-12*h_tot1(i)) then
735  write(mesg,'(ES11.4," became ",ES11.4," diff ",ES11.4)') &
736  h_tot1(i), h_tot2(i), (h_tot1(i) - h_tot2(i))
737  call mom_error(warning, "regularize_surface: Mass non-conservation."//&
738  trim(mesg), .true.)
739  fatal_error = .true.
740  endif
741  if (abs(th_tot1(i) - th_tot2(i)) > 1e-12*(th_tot1(i)+10.0*h_tot1(i))) then
742  write(mesg,'(ES11.4," became ",ES11.4," diff ",ES11.4," int diff ",ES11.4)') &
743  th_tot1(i), th_tot2(i), (th_tot1(i) - th_tot2(i)), (th_tot1(i) - th_tot3(i))
744  call mom_error(warning, "regularize_surface: Heat non-conservation."//&
745  trim(mesg), .true.)
746  fatal_error = .true.
747  endif
748  if (abs(sh_tot1(i) - sh_tot2(i)) > 1e-12*(sh_tot1(i)+10.0*h_tot1(i))) then
749  write(mesg,'(ES11.4," became ",ES11.4," diff ",ES11.4," int diff ",ES11.4)') &
750  sh_tot1(i), sh_tot2(i), (sh_tot1(i) - sh_tot2(i)), (sh_tot1(i) - sh_tot3(i))
751  call mom_error(warning, "regularize_surface: Salinity non-conservation."//&
752  trim(mesg), .true.)
753  fatal_error = .true.
754  endif
755  if (fatal_error) then
756  write(mesg,'("Error at lat/lon ",2(ES11.4))') g%geoLatT(i,j), g%geoLonT(i,j)
757  call mom_error(fatal, "regularize_surface: Terminating with fatal error. "//&
758  trim(mesg))
759  endif
760  endif ; enddo
761  endif
762 
763  endif ; enddo ! j-loop.
764 
765  if (cs%id_def_rat > 0) call post_data(cs%id_def_rat, def_rat_h, cs%diag)
766 
767 #ifdef DEBUG_CODE
768  if (cs%id_e1 > 0) call post_data(cs%id_e1, e, cs%diag)
769  if (cs%id_def_rat_u > 0) call post_data(cs%id_def_rat_u, def_rat_u, cs%diag)
770  if (cs%id_def_rat_u_1b > 0) call post_data(cs%id_def_rat_u_1b, def_rat_u_1b, cs%diag)
771  if (cs%id_def_rat_v > 0) call post_data(cs%id_def_rat_v, def_rat_v, cs%diag)
772  if (cs%id_def_rat_v_1b > 0) call post_data(cs%id_def_rat_v_1b, def_rat_v_1b, cs%diag)
773 
774  if ((cs%id_def_rat_2 > 0) .or. &
775  (cs%id_def_rat_u_2 > 0) .or. (cs%id_def_rat_u_2b > 0) .or. &
776  (cs%id_def_rat_v_2 > 0) .or. (cs%id_def_rat_v_2b > 0) ) then
777  do j=js-1,je+1 ; do i=is-1,ie+1
778  e(i,j,1) = 0.0
779  enddo ; enddo
780  do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
781  e(i,j,k+1) = e(i,j,k) - h(i,j,k)
782  enddo ; enddo ; enddo
783 
784  call find_deficit_ratios(e, def_rat_u_2, def_rat_v_2, g, gv, cs, def_rat_u_2b, def_rat_v_2b, h=h)
785 
786  ! Determine which columns are problematic
787  do j=js,je ; do i=is,ie
788  def_rat_h2(i,j) = max(def_rat_u_2(i-1,j), def_rat_u_2(i,j), &
789  def_rat_v_2(i,j-1), def_rat_v_2(i,j))
790  enddo ; enddo
791 
792  if (cs%id_def_rat_2 > 0) call post_data(cs%id_def_rat_2, def_rat_h2, cs%diag)
793  if (cs%id_e2 > 0) call post_data(cs%id_e2, e, cs%diag)
794  if (cs%id_def_rat_u_2 > 0) call post_data(cs%id_def_rat_u_2, def_rat_u_2, cs%diag)
795  if (cs%id_def_rat_u_2b > 0) call post_data(cs%id_def_rat_u_2b, def_rat_u_2b, cs%diag)
796  if (cs%id_def_rat_v_2 > 0) call post_data(cs%id_def_rat_v_2, def_rat_v_2, cs%diag)
797  if (cs%id_def_rat_v_2b > 0) call post_data(cs%id_def_rat_v_2b, def_rat_v_2b, cs%diag)
798  endif
799 #endif
800 
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:19
Here is the call graph for this function:
Here is the caller graph for this function:

Variable Documentation

◆ id_clock_eos

integer mom_regularize_layers::id_clock_eos
private

Definition at line 111 of file MOM_regularize_layers.F90.

Referenced by regularize_layers_init(), and regularize_surface().

◆ id_clock_pass

integer mom_regularize_layers::id_clock_pass

Definition at line 111 of file MOM_regularize_layers.F90.

Referenced by regularize_layers(), and regularize_layers_init().

111 integer :: id_clock_pass, id_clock_eos