user_change_diffusivity.F90

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