MOM6
user_change_diffusivity.F90
Go to the documentation of this file.
1 !> Increments the diapycnal diffusivity in a specified band of latitudes and densities.
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
6 use mom_diag_mediator, only : diag_ctrl, time_type
7 use mom_error_handler, only : mom_error, is_root_pe, fatal, warning, note
9 use mom_grid, only : ocean_grid_type
13 use mom_eos, only : calculate_density
14 
15 implicit none ; private
16 
17 #include <MOM_memory.h>
18 
20 
21 ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
22 ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
23 ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
24 ! vary with the Boussinesq approximation, the Boussinesq variant is given first.
25 
26 !> Control structure for user_change_diffusivity
27 type, public :: user_change_diff_cs ; private
28  real :: kd_add !< The scale of a diffusivity that is added everywhere
29  !! without any filtering or scaling [Z2 T-1 ~> m2 s-1].
30  real :: lat_range(4) !< 4 values that define the latitude range over which
31  !! a diffusivity scaled by Kd_add is added [degLat].
32  real :: rho_range(4) !< 4 values that define the coordinate potential
33  !! density range over which a diffusivity scaled by
34  !! Kd_add is added [kg m-3].
35  logical :: use_abs_lat !< If true, use the absolute value of latitude when
36  !! setting lat_range.
37  type(diag_ctrl), pointer :: diag => null() !< A structure that is used to
38  !! regulate the timing of diagnostic output.
39 end type user_change_diff_cs
40 
41 contains
42 
43 !> This subroutine provides an interface for a user to use to modify the
44 !! main code to alter the diffusivities as needed. The specific example
45 !! implemented here augments the diffusivity for a specified range of latitude
46 !! and coordinate potential density.
47 subroutine user_change_diff(h, tv, G, GV, CS, Kd_lay, Kd_int, T_f, S_f, Kd_int_add)
48  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
49  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
50  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2].
51  type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers
52  !! to any available thermodynamic
53  !! fields. Absent fields have NULL ptrs.
54  type(user_change_diff_cs), pointer :: cs !< This module's control structure.
55  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), optional, intent(inout) :: kd_lay !< The diapycnal diffusivity of
56  !! each layer [Z2 T-1 ~> m2 s-1].
57  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), optional, intent(inout) :: kd_int !< The diapycnal diffusivity
58  !! at each interface [Z2 T-1 ~> m2 s-1].
59  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), optional, intent(in) :: t_f !< Temperature with massless
60  !! layers filled in vertically [degC].
61  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), optional, intent(in) :: s_f !< Salinity with massless
62  !! layers filled in vertically [ppt].
63  real, dimension(:,:,:), optional, pointer :: kd_int_add !< The diapycnal
64  !! diffusivity that is being added at
65  !! each interface [Z2 T-1 ~> m2 s-1].
66  ! Local variables
67  real :: rcv(szi_(g),szk_(g)) ! The coordinate density in layers [kg m-3].
68  real :: p_ref(szi_(g)) ! An array of tv%P_Ref pressures.
69  real :: rho_fn ! The density dependence of the input function, 0-1 [nondim].
70  real :: lat_fn ! The latitude dependence of the input function, 0-1 [nondim].
71  logical :: use_eos ! If true, density is calculated from T & S using an
72  ! equation of state.
73  logical :: store_kd_add ! Save the added diffusivity as a diagnostic if true.
74  integer :: i, j, k, is, ie, js, je, nz
75  integer :: isd, ied, jsd, jed
76 
77  real :: kappa_fill ! diffusivity used to fill massless layers
78  real :: dt_fill ! timestep used to fill massless layers
79  character(len=200) :: mesg
80 
81  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
82  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
83 
84  if (.not.associated(cs)) call mom_error(fatal,"user_set_diffusivity: "//&
85  "Module must be initialized before it is used.")
86 
87  use_eos = associated(tv%eqn_of_state)
88  if (.not.use_eos) return
89  store_kd_add = .false.
90  if (present(kd_int_add)) store_kd_add = associated(kd_int_add)
91 
92  if (.not.range_ok(cs%lat_range)) then
93  write(mesg, '(4(1pe15.6))') cs%lat_range(1:4)
94  call mom_error(fatal, "user_set_diffusivity: bad latitude range: \n "//&
95  trim(mesg))
96  endif
97  if (.not.range_ok(cs%rho_range)) then
98  write(mesg, '(4(1pe15.6))') cs%rho_range(1:4)
99  call mom_error(fatal, "user_set_diffusivity: bad density range: \n "//&
100  trim(mesg))
101  endif
102 
103  if (store_kd_add) kd_int_add(:,:,:) = 0.0
104 
105  do i=is,ie ; p_ref(i) = tv%P_Ref ; enddo
106  do j=js,je
107  if (present(t_f) .and. present(s_f)) then
108  do k=1,nz
109  call calculate_density(t_f(:,j,k),s_f(:,j,k),p_ref,rcv(:,k),&
110  is,ie-is+1,tv%eqn_of_state)
111  enddo
112  else
113  do k=1,nz
114  call calculate_density(tv%T(:,j,k),tv%S(:,j,k),p_ref,rcv(:,k),&
115  is,ie-is+1,tv%eqn_of_state)
116  enddo
117  endif
118 
119  if (present(kd_lay)) then
120  do k=1,nz ; do i=is,ie
121  if (cs%use_abs_lat) then
122  lat_fn = val_weights(abs(g%geoLatT(i,j)), cs%lat_range)
123  else
124  lat_fn = val_weights(g%geoLatT(i,j), cs%lat_range)
125  endif
126  rho_fn = val_weights(rcv(i,k), cs%rho_range)
127  if (rho_fn * lat_fn > 0.0) &
128  kd_lay(i,j,k) = kd_lay(i,j,k) + cs%Kd_add * rho_fn * lat_fn
129  enddo ; enddo
130  endif
131  if (present(kd_int)) then
132  do k=2,nz ; do i=is,ie
133  if (cs%use_abs_lat) then
134  lat_fn = val_weights(abs(g%geoLatT(i,j)), cs%lat_range)
135  else
136  lat_fn = val_weights(g%geoLatT(i,j), cs%lat_range)
137  endif
138  ! rho_int = 0.5*(Rcv(i,k-1) + Rcv(i,k))
139  rho_fn = val_weights( 0.5*(rcv(i,k-1) + rcv(i,k)), cs%rho_range)
140  if (rho_fn * lat_fn > 0.0) then
141  kd_int(i,j,k) = kd_int(i,j,k) + cs%Kd_add * rho_fn * lat_fn
142  if (store_kd_add) kd_int_add(i,j,k) = cs%Kd_add * rho_fn * lat_fn
143  endif
144  enddo ; enddo
145  endif
146  enddo
147 
148 end subroutine user_change_diff
149 
150 !> This subroutine checks whether the 4 values of range are in ascending order.
151 function range_ok(range) result(OK)
152  real, dimension(4), intent(in) :: range !< Four values to check.
153  logical :: ok !< Return value.
154 
155  ok = ((range(1) <= range(2)) .and. (range(2) <= range(3)) .and. &
156  (range(3) <= range(4)))
157 
158 end function range_ok
159 
160 !> This subroutine returns a value that goes smoothly from 0 to 1, stays
161 !! at 1, and then goes smoothly back to 0 at the four values of range. The
162 !! transitions are cubic, and have zero first derivatives where the curves
163 !! hit 0 and 1. The values in range must be in ascending order, as can be
164 !! checked by calling range_OK.
165 function val_weights(val, range) result(ans)
166  real, intent(in) :: val !< Value for which we need an answer.
167  real, dimension(4), intent(in) :: range !< Range over which the answer is non-zero.
168  real :: ans !< Return value.
169  ! Local variables
170  real :: x ! A nondimensional number between 0 and 1.
171 
172  ans = 0.0
173  if ((val > range(1)) .and. (val < range(4))) then
174  if (val < range(2)) then
175  ! x goes from 0 to 1; ans goes from 0 to 1, with 0 derivatives at the ends.
176  x = (val - range(1)) / (range(2) - range(1))
177  ans = x**2 * (3.0 - 2.0 * x)
178  elseif (val > range(3)) then
179  ! x goes from 0 to 1; ans goes from 0 to 1, with 0 derivatives at the ends.
180  x = (range(4) - val) / (range(4) - range(3))
181  ans = x**2 * (3.0 - 2.0 * x)
182  else
183  ans = 1.0
184  endif
185  endif
186 
187 end function val_weights
188 
189 !> Set up the module control structure.
190 subroutine user_change_diff_init(Time, G, GV, US, param_file, diag, CS)
191  type(time_type), intent(in) :: time !< The current model time.
192  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
193  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
194  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
195  type(param_file_type), intent(in) :: param_file !< A structure indicating the
196  !! open file to parse for
197  !! model parameter values.
198  type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to
199  !! regulate diagnostic output.
200  type(user_change_diff_cs), pointer :: cs !< A pointer that is set to
201  !! point to the control
202  !! structure for this module.
203 
204 ! This include declares and sets the variable "version".
205 #include "version_variable.h"
206  character(len=40) :: mdl = "user_set_diffusivity" ! This module's name.
207  character(len=200) :: mesg
208  integer :: i, j, is, ie, js, je
209 
210  if (associated(cs)) then
211  call mom_error(warning, "diabatic_entrain_init called with an associated "// &
212  "control structure.")
213  return
214  endif
215  allocate(cs)
216 
217  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
218 
219  cs%diag => diag
220 
221  ! Read all relevant parameters and write them to the model log.
222  call log_version(param_file, mdl, version, "")
223  call get_param(param_file, mdl, "USER_KD_ADD", cs%Kd_add, &
224  "A user-specified additional diffusivity over a range of "//&
225  "latitude and density.", default=0.0, units="m2 s-1", &
226  scale=us%m2_s_to_Z2_T)
227  if (cs%Kd_add /= 0.0) then
228  call get_param(param_file, mdl, "USER_KD_ADD_LAT_RANGE", cs%lat_range(:), &
229  "Four successive values that define a range of latitudes "//&
230  "over which the user-specified extra diffusivity is "//&
231  "applied. The four values specify the latitudes at "//&
232  "which the extra diffusivity starts to increase from 0, "//&
233  "hits its full value, starts to decrease again, and is "//&
234  "back to 0.", units="degree", default=-1.0e9)
235  call get_param(param_file, mdl, "USER_KD_ADD_RHO_RANGE", cs%rho_range(:), &
236  "Four successive values that define a range of potential "//&
237  "densities over which the user-given extra diffusivity "//&
238  "is applied. The four values specify the density at "//&
239  "which the extra diffusivity starts to increase from 0, "//&
240  "hits its full value, starts to decrease again, and is "//&
241  "back to 0.", units="kg m-3", default=-1.0e9)
242  call get_param(param_file, mdl, "USER_KD_ADD_USE_ABS_LAT", cs%use_abs_lat, &
243  "If true, use the absolute value of latitude when "//&
244  "checking whether a point fits into range of latitudes.", &
245  default=.false.)
246  endif
247 
248  if (.not.range_ok(cs%lat_range)) then
249  write(mesg, '(4(1pe15.6))') cs%lat_range(1:4)
250  call mom_error(fatal, "user_set_diffusivity: bad latitude range: \n "//&
251  trim(mesg))
252  endif
253  if (.not.range_ok(cs%rho_range)) then
254  write(mesg, '(4(1pe15.6))') cs%rho_range(1:4)
255  call mom_error(fatal, "user_set_diffusivity: bad density range: \n "//&
256  trim(mesg))
257  endif
258 
259 end subroutine user_change_diff_init
260 
261 !> Clean up the module control structure.
262 subroutine user_change_diff_end(CS)
263  type(user_change_diff_cs), pointer :: cs !< A pointer that is set to point to the control
264  !! structure for this module.
265 
266  if (associated(cs)) deallocate(cs)
267 
268 end subroutine user_change_diff_end
269 
270 end module user_change_diffusivity
mom_verticalgrid
Provides a transparent vertical ocean grid type and supporting routines.
Definition: MOM_verticalGrid.F90:2
mom_file_parser::log_version
An overloaded interface to log version information about modules.
Definition: MOM_file_parser.F90:109
mom_diag_mediator
The subroutines here provide convenient wrappers to the fms diag_manager interfaces with additional d...
Definition: MOM_diag_mediator.F90:3
mom_variables::thermo_var_ptrs
Pointers to an assortment of thermodynamic fields that may be available, including potential temperat...
Definition: MOM_variables.F90:78
mom_eos
Provides subroutines for quantities specific to the equation of state.
Definition: MOM_EOS.F90:2
user_change_diffusivity::user_change_diff_end
subroutine, public user_change_diff_end(CS)
Clean up the module control structure.
Definition: user_change_diffusivity.F90:263
mom_file_parser::param_file_type
A structure that can be parsed to read and document run-time parameters.
Definition: MOM_file_parser.F90:54
mom_file_parser::get_param
An overloaded interface to read and log the values of various types of parameters.
Definition: MOM_file_parser.F90:102
user_change_diffusivity::user_change_diff
subroutine, public user_change_diff(h, tv, G, GV, CS, Kd_lay, Kd_int, T_f, S_f, Kd_int_add)
This subroutine provides an interface for a user to use to modify the main code to alter the diffusiv...
Definition: user_change_diffusivity.F90:48
mom_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
user_change_diffusivity::user_change_diff_init
subroutine, public user_change_diff_init(Time, G, GV, US, param_file, diag, CS)
Set up the module control structure.
Definition: user_change_diffusivity.F90:191
mom_variables::vertvisc_type
Vertical viscosities, drag coefficients, and related fields.
Definition: MOM_variables.F90:196
mom_verticalgrid::verticalgrid_type
Describes the vertical ocean grid, including unit conversion factors.
Definition: MOM_verticalGrid.F90:24
mom_variables
Provides transparent structures with groups of MOM6 variables and supporting routines.
Definition: MOM_variables.F90:2
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
user_change_diffusivity
Increments the diapycnal diffusivity in a specified band of latitudes and densities.
Definition: user_change_diffusivity.F90:2
user_change_diffusivity::user_change_diff_cs
Control structure for user_change_diffusivity.
Definition: user_change_diffusivity.F90:27
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.F90:2
mom_unit_scaling
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Definition: MOM_unit_scaling.F90:2
mom_error_handler::is_root_pe
logical function, public is_root_pe()
This returns .true. if the current PE is the root PE.
Definition: MOM_error_handler.F90:44
user_change_diffusivity::range_ok
logical function range_ok(range)
This subroutine checks whether the 4 values of range are in ascending order.
Definition: user_change_diffusivity.F90:152
mom_variables::p3d
A structure for creating arrays of pointers to 3D arrays.
Definition: MOM_variables.F90:28
mom_error_handler::mom_error
subroutine, public mom_error(level, message, all_print)
This provides a convenient interface for writing an mpp_error message with run-time filter based on a...
Definition: MOM_error_handler.F90:72
user_change_diffusivity::val_weights
real function val_weights(val, range)
This subroutine returns a value that goes smoothly from 0 to 1, stays at 1, and then goes smoothly ba...
Definition: user_change_diffusivity.F90:166
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:26
mom_diag_mediator::diag_ctrl
The following data type a list of diagnostic fields an their variants, as well as variables that cont...
Definition: MOM_diag_mediator.F90:239
mom_eos::calculate_density
Calculates density of sea water from T, S and P.
Definition: MOM_EOS.F90:60