-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathuser_change_diffusivity.F90
328 lines (285 loc) · 15.9 KB
/
user_change_diffusivity.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
!> Increments the diapycnal diffusivity in a specified band of latitudes and densities.
module user_change_diffusivity
! This file is part of MOM6. See LICENSE.md for the license.
use MOM_diag_mediator, only : diag_ctrl, time_type
use MOM_error_handler, only : MOM_error, is_root_pe, FATAL, WARNING, NOTE
use MOM_file_parser, only : get_param, log_version, param_file_type
use MOM_grid, only : ocean_grid_type
use MOM_unit_scaling, only : unit_scale_type
use MOM_variables, only : thermo_var_ptrs, vertvisc_type, p3d
use MOM_verticalGrid, only : verticalGrid_type
use MOM_EOS, only : calculate_density, EOS_domain
use MOM_forcing_type, only : forcing
implicit none ; private
#include <MOM_memory.h>
public user_change_diff, user_change_diff_init, user_change_diff_end
! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
! vary with the Boussinesq approximation, the Boussinesq variant is given first.
!> Control structure for user_change_diffusivity
type, public :: user_change_diff_CS ; private
logical :: initialized = .false. !< True if this control structure has been initialized.
real :: Kd_add !< The scale of a diffusivity that is added everywhere
!! without any filtering or scaling [Z2 T-1 ~> m2 s-1].
real :: lat_range(4) !< 4 values that define the latitude range over which
!! a diffusivity scaled by Kd_add is added [degrees_N].
real :: rho_range(4) !< 4 values that define the coordinate potential
!! density range over which a diffusivity scaled by
!! Kd_add is added [R ~> kg m-3].
real :: dep_range(4) !< 4 values that define the coordinate depth range [L ~> m]
!! over which a diffusivity scaled by
!! Kd_add is added [R ~> kg m-3].
logical :: use_abs_lat !< If true, use the absolute value of latitude when
!! setting lat_range.
logical :: use_dep_range !< If true, use the actual value of depth when
!! setting dep_range.
type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to
!! regulate the timing of diagnostic output.
end type user_change_diff_CS
contains
!> This subroutine provides an interface for a user to use to modify the
!! main code to alter the diffusivities as needed. The specific example
!! implemented here augments the diffusivity for a specified range of latitude
!! and coordinate potential density.
subroutine user_change_diff(fluxes, h, tv, G, GV, US, CS, Kd_lay, Kd_int, T_f, S_f, Kd_int_add)
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2].
type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers
!! to any available thermodynamic
!! fields. Absent fields have NULL ptrs.
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
type(user_change_diff_CS), pointer :: CS !< This module's control structure.
type(forcing), intent(in) :: fluxes !< Surface fluxes container
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), optional, intent(inout) :: Kd_lay !< The diapycnal diffusivity of
!! each layer [Z2 T-1 ~> m2 s-1].
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), optional, intent(inout) :: Kd_int !< The diapycnal diffusivity
!! at each interface [Z2 T-1 ~> m2 s-1].
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), optional, intent(in) :: T_f !< Temperature with massless
!! layers filled in vertically [C ~> degC].
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), optional, intent(in) :: S_f !< Salinity with massless
!! layers filled in vertically [S ~> ppt].
real, dimension(:,:,:), optional, pointer :: Kd_int_add !< The diapycnal
!! diffusivity that is being added at
!! each interface [Z2 T-1 ~> m2 s-1].
! Local variables
real :: Rcv(SZI_(G),SZK_(GV)) ! The coordinate density in layers [R ~> kg m-3].
real :: p_ref(SZI_(G)) ! An array of tv%P_Ref pressures [R L2 T-2 ~> Pa].
real :: rho_fn ! The density dependence of the input function, 0-1 [nondim].
real :: dep_fn ! The depth dependence of the input function, 0-1 [nondim].
real :: lat_fn ! The latitude dependence of the input function, 0-1 [nondim].
logical :: use_EOS ! If true, density is calculated from T & S using an
! equation of state.
logical :: store_Kd_add ! Save the added diffusivity as a diagnostic if true.
integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
integer :: i, j, k, is, ie, js, je, nz
integer :: isd, ied, jsd, jed
real :: zmom(SZI_(G),SZJ_(G),SZK_(GV)) ! An array of vertical grid depth [L ~> m]
character(len=200) :: mesg
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
!!Vertical Grid Depth (zmom) using layer thickness value (h)
do k=1,nz
if ( k .eq. 1 ) then
zmom(:,:,k) = 0.5*h(:,:,k)
else
zmom(:,:,k) = zmom(:,:,k-1) + 0.5*(h(:,:,k-1)+h(:,:,k))
endif
enddo
if (.not.associated(CS)) call MOM_error(FATAL,"user_set_diffusivity: "//&
"Module must be initialized before it is used.")
if (.not.CS%initialized) call MOM_error(FATAL,"user_set_diffusivity: "//&
"Module must be initialized before it is used.")
use_EOS = associated(tv%eqn_of_state)
if (.not.use_EOS) return
store_Kd_add = .false.
if (present(Kd_int_add)) store_Kd_add = associated(Kd_int_add)
if (.not.range_OK(CS%lat_range)) then
write(mesg, '(4(1pe15.6))') CS%lat_range(1:4)
call MOM_error(FATAL, "user_set_diffusivity: bad latitude range: \n "//&
trim(mesg))
endif
if (.not.range_OK(CS%rho_range)) then
write(mesg, '(4(1pe15.6))') CS%rho_range(1:4)
call MOM_error(FATAL, "user_set_diffusivity: bad density range: \n "//&
trim(mesg))
endif
if (.not.range_OK(CS%dep_range)) then
write(mesg, '(4(1pe15.6))') CS%dep_range(1:4)
call MOM_error(FATAL, "user_set_diffusivity: bad depth range: \n "//&
trim(mesg))
endif
if (store_Kd_add) Kd_int_add(:,:,:) = 0.0
do i=is,ie ; p_ref(i) = tv%P_Ref ; enddo
EOSdom(:) = EOS_domain(G%HI)
do j=js,je !! J loop starts
if (present(T_f) .and. present(S_f)) then
do k=1,nz
call calculate_density(T_f(:,j,k), S_f(:,j,k), p_ref, Rcv(:,k), tv%eqn_of_state, EOSdom)
enddo
else
do k=1,nz
call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_ref, Rcv(:,k), tv%eqn_of_state, EOSdom)
enddo
endif
if (present(Kd_lay)) then
do k=1,nz ; do i=is,ie
if (CS%use_abs_lat) then
lat_fn = val_weights(abs(G%geoLatT(i,j)), CS%lat_range)
else
lat_fn = val_weights(G%geoLatT(i,j), CS%lat_range)
endif
rho_fn = val_weights(Rcv(i,k), CS%rho_range)
if ((CS%use_dep_range) .and. (associated(fluxes%lrunoff))) then
if ((zmom(i,j,k) .lt. maxval(CS%dep_range)) .and. (fluxes%lrunoff(i,j) .ne. 0.0)) then
Kd_lay(i,j,k) = max(Kd_lay(i,j,k),CS%Kd_add)
endif
else
if (rho_fn * lat_fn > 0.0) &
Kd_lay(i,j,k) = Kd_lay(i,j,k) + CS%Kd_add * rho_fn * lat_fn
endif
enddo ; enddo
endif
if (present(Kd_int)) then
do K=2,nz ; do i=is,ie
if (CS%use_abs_lat) then
lat_fn = val_weights(abs(G%geoLatT(i,j)), CS%lat_range)
else
lat_fn = val_weights(G%geoLatT(i,j), CS%lat_range)
endif
rho_fn = val_weights( 0.5*(Rcv(i,k-1) + Rcv(i,k)), CS%rho_range)
dep_fn = val_weights(zmom(i,j,k-1),CS%dep_range)
if ((CS%use_dep_range) .and. (associated(fluxes%lrunoff))) then
if ((zmom(i,j,k) .lt. maxval(CS%dep_range)) .and. (fluxes%lrunoff(i,j) .ne. 0.0)) then
Kd_int(i,j,k) = max(Kd_int(i,j,k),CS%Kd_add)
if (store_Kd_add) Kd_int_add(i,j,K) = CS%Kd_add * dep_fn * lat_fn
endif
else
if (rho_fn * lat_fn > 0.0) then
Kd_int(i,j,K) = Kd_int(i,j,K) + CS%Kd_add * rho_fn * lat_fn
if (store_Kd_add) Kd_int_add(i,j,K) = CS%Kd_add * rho_fn * lat_fn
endif
endif
enddo ; enddo
endif
enddo !! J loop ends
end subroutine user_change_diff
!> This subroutine checks whether the 4 values of range are in ascending order.
function range_OK(range) result(OK)
real, dimension(4), intent(in) :: range !< Four values to check [arbitrary]
logical :: OK !< Return value.
OK = ((range(1) <= range(2)) .and. (range(2) <= range(3)) .and. &
(range(3) <= range(4)))
end function range_OK
!> This subroutine returns a value that goes smoothly from 0 to 1, stays
!! at 1, and then goes smoothly back to 0 at the four values of range. The
!! transitions are cubic, and have zero first derivatives where the curves
!! hit 0 and 1. The values in range must be in ascending order, as can be
!! checked by calling range_OK.
function val_weights(val, range) result(ans)
real, intent(in) :: val !< Value for which we need an answer [arbitrary units].
real, dimension(4), intent(in) :: range !< Range over which the answer is non-zero [arbitrary units].
real :: ans !< Return value [nondim].
! Local variables
real :: x ! A nondimensional number between 0 and 1 [nondim].
ans = 0.0
if ((val > range(1)) .and. (val < range(4))) then
if (val < range(2)) then
! x goes from 0 to 1; ans goes from 0 to 1, with 0 derivatives at the ends.
x = (val - range(1)) / (range(2) - range(1))
ans = x**2 * (3.0 - 2.0 * x)
elseif (val > range(3)) then
! x goes from 0 to 1; ans goes from 0 to 1, with 0 derivatives at the ends.
x = (range(4) - val) / (range(4) - range(3))
ans = x**2 * (3.0 - 2.0 * x)
else
ans = 1.0
endif
endif
end function val_weights
!> Set up the module control structure.
subroutine user_change_diff_init(Time, G, GV, US, param_file, diag, CS)
type(time_type), intent(in) :: Time !< The current model time.
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
type(param_file_type), intent(in) :: param_file !< A structure indicating the
!! open file to parse for
!! model parameter values.
type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to
!! regulate diagnostic output.
type(user_change_diff_CS), pointer :: CS !< A pointer that is set to
!! point to the control
!! structure for this module.
! This include declares and sets the variable "version".
# include "version_variable.h"
character(len=40) :: mdl = "user_set_diffusivity" ! This module's name.
character(len=200) :: mesg
if (associated(CS)) then
call MOM_error(WARNING, "diabatic_entrain_init called with an associated "// &
"control structure.")
return
endif
allocate(CS)
CS%initialized = .true.
CS%diag => diag
! Read all relevant parameters and write them to the model log.
call log_version(param_file, mdl, version, "")
call get_param(param_file, mdl, "USER_KD_ADD", CS%Kd_add, &
"A user-specified additional diffusivity over a range of "//&
"latitude and density.", default=0.0, units="m2 s-1", scale=US%m2_s_to_Z2_T)
if (CS%Kd_add /= 0.0) then
call get_param(param_file, mdl, "USER_KD_ADD_LAT_RANGE", CS%lat_range(:), &
"Four successive values that define a range of latitudes "//&
"over which the user-specified extra diffusivity is "//&
"applied. The four values specify the latitudes at "//&
"which the extra diffusivity starts to increase from 0, "//&
"hits its full value, starts to decrease again, and is "//&
"back to 0.", units="degrees_N", default=-1.0e9)
call get_param(param_file, mdl, "USER_KD_ADD_RHO_RANGE", CS%rho_range(:), &
"Four successive values that define a range of potential "//&
"densities over which the user-given extra diffusivity "//&
"is applied. The four values specify the density at "//&
"which the extra diffusivity starts to increase from 0, "//&
"hits its full value, starts to decrease again, and is "//&
"back to 0.", units="kg m-3", default=-1.0e9, scale=US%kg_m3_to_R)
call get_param(param_file, mdl, "USER_KD_ADD_DEP_RANGE", CS%dep_range(:), &
"Four successive values that define a range of depth "//&
"over which the user-given extra diffusivity "//&
"is applied. The four values specify the depth at "//&
"which the extra diffusivity starts to increase from 0, "//&
"hits its full value, starts to decrease again, and is "//&
"back to 0.", units="m", default=-1.0e9, scale=US%m_to_Z)
call get_param(param_file, mdl, "USER_KD_ADD_USE_ABS_LAT", CS%use_abs_lat, &
"If true, use the absolute value of latitude when "//&
"checking whether a point fits into range of latitudes.", &
default=.false.)
call get_param(param_file, mdl, "USER_KD_ADD_USE_DEPTH_RANGE", CS%use_dep_range, &
"If true, use the actual value of depth when "//&
"checking whether a point fits into range of depths.", &
default=.false.)
endif
if (.not.range_OK(CS%lat_range)) then
write(mesg, '(4(1pe15.6))') CS%lat_range(1:4)
call MOM_error(FATAL, "user_set_diffusivity: bad latitude range: \n "//&
trim(mesg))
endif
if (.not.range_OK(CS%rho_range)) then
write(mesg, '(4(1pe15.6))') CS%rho_range(1:4)
call MOM_error(FATAL, "user_set_diffusivity: bad density range: \n "//&
trim(mesg))
endif
if (.not.range_OK(CS%dep_range)) then
write(mesg, '(4(1pe15.6))') CS%dep_range(1:4)
call MOM_error(FATAL, "user_set_diffusivity: bad depth range: \n "//&
trim(mesg))
endif
end subroutine user_change_diff_init
!> Clean up the module control structure.
subroutine user_change_diff_end(CS)
type(user_change_diff_CS), pointer :: CS !< A pointer that is set to point to the control
!! structure for this module.
if (associated(CS)) deallocate(CS)
end subroutine user_change_diff_end
end module user_change_diffusivity