dyed_channel_initialization.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!> Initialization for the dyed_channel configuration
7
9use mom_error_handler, only : mom_mesg, mom_error, fatal, warning, is_root_pe
10use mom_file_parser, only : get_param, log_version, param_file_type
11use mom_get_input, only : directories
12use mom_grid, only : ocean_grid_type
13use mom_open_boundary, only : ocean_obc_type, obc_none
14use mom_open_boundary, only : obc_direction_w, obc_direction_n, obc_direction_s, obc_direction_e
17use mom_time_manager, only : time_type, time_to_real
18use mom_tracer_registry, only : tracer_registry_type, tracer_name_lookup
19use mom_tracer_registry, only : tracer_type
23
24implicit none ; private
25
26#include <MOM_memory.h>
27
30
31!> Control structure for dyed-channel open boundaries.
32type, public :: dyed_channel_obc_cs ; private
33 real :: zonal_flow = 8.57 !< Mean inflow [L T-1 ~> m s-1]
34 real :: tidal_amp = 0.0 !< Sloshing amplitude [L T-1 ~> m s-1]
35 real :: frequency = 0.0 !< Sloshing frequency [T-1 ~> s-1]
36 logical :: obc_transport_bug !< If true and specified open boundary conditions are being
37 !! used, use a 1 m (if Boussienesq) or 1 kg m-2 layer thickness
38 !! instead of the actual thickness.
40
41integer :: ntr = 0 !< Number of dye tracers
42 !! \todo This is a module variable. Move this variable into the control structure.
43
44contains
45
46!> Add dyed channel to OBC registry.
47logical function register_dyed_channel_obc(param_file, CS, US, OBC_Reg)
48 type(param_file_type), intent(in) :: param_file !< parameter file.
49 type(dyed_channel_obc_cs), pointer :: cs !< Dyed channel control structure.
50 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
51 type(obc_registry_type), pointer :: obc_reg !< OBC registry.
52
53 ! Local variables
54 logical :: enable_bugs ! If true, the defaults for recently added bug-fix flags are set to
55 ! recreate the bugs, or if false bugs are only used if actively selected.
56 character(len=32) :: casename = "dyed channel" ! This case's name.
57 character(len=40) :: mdl = "register_dyed_channel_OBC" ! This subroutine's name.
58
59 if (associated(cs)) then
60 call mom_error(warning, "register_dyed_channel_OBC called with an "// &
61 "associated control structure.")
62 return
63 endif
64 allocate(cs)
65
66 call get_param(param_file, mdl, "CHANNEL_MEAN_FLOW", cs%zonal_flow, &
67 "Mean zonal flow imposed at upstream open boundary.", &
68 units="m/s", default=8.57, scale=us%m_s_to_L_T)
69 call get_param(param_file, mdl, "CHANNEL_TIDAL_AMP", cs%tidal_amp, &
70 "Sloshing amplitude imposed at upstream open boundary.", &
71 units="m/s", default=0.0, scale=us%m_s_to_L_T)
72 call get_param(param_file, mdl, "CHANNEL_FLOW_FREQUENCY", cs%frequency, &
73 "Frequency of oscillating zonal flow.", &
74 units="s-1", default=0.0, scale=us%T_to_s)
75 call get_param(param_file, mdl, "ENABLE_BUGS_BY_DEFAULT", enable_bugs, &
76 default=.true., do_not_log=.true.) ! This is logged from MOM.F90.
77 call get_param(param_file, mdl, "CHANNEL_FLOW_OBC_TRANSPORT_BUG", cs%OBC_transport_bug, &
78 "If true and specified open boundary conditions are being used, use a 1 m "//&
79 "(if Boussienesq) or 1 kg m-2 layer thickness instead of the actual thickness.", &
80 default=enable_bugs)
81
82 ! Register the open boundaries.
83 call register_obc(casename, param_file, obc_reg)
85
87
88!> Clean up the dyed_channel OBC from registry.
89subroutine dyed_channel_obc_end(CS)
90 type(dyed_channel_obc_cs), pointer :: cs !< Dyed channel control structure.
91
92 if (associated(cs)) then
93 deallocate(cs)
94 endif
95end subroutine dyed_channel_obc_end
96
97!> This subroutine sets the dye and flow properties at open boundary conditions.
98subroutine dyed_channel_set_obc_tracer_data(OBC, G, GV, param_file, tr_Reg)
99 type(ocean_obc_type), pointer :: obc !< This open boundary condition type specifies
100 !! whether, where, and what open boundary
101 !! conditions are used.
102 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
103 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
104 type(param_file_type), intent(in) :: param_file !< A structure indicating the open file
105 !! to parse for model parameter values.
106 type(tracer_registry_type), pointer :: tr_reg !< Tracer registry.
107 ! Local variables
108 character(len=40) :: mdl = "dyed_channel_set_OBC_tracer_data" ! This subroutine's name.
109 character(len=80) :: name, longname
110 integer :: m, n, ntr_id
111 real :: dye ! Inflow dye concentrations [arbitrary]
112 type(tracer_type), pointer :: tr_ptr => null()
113
114 if (.not.associated(obc)) call mom_error(fatal, 'dyed_channel_initialization.F90: '// &
115 'dyed_channel_set_OBC_data() was called but OBC type was not initialized!')
116
117 call get_param(param_file, mdl, "NUM_DYE_TRACERS", ntr, &
118 "The number of dye tracers in this run. Each tracer "//&
119 "should have a separate boundary segment.", default=0, &
120 do_not_log=.true.)
121
122 if (obc%number_of_segments < ntr) then
123 call mom_error(warning, "Error in dyed_obc segment setup")
124 return !!! Need a better error message here
125 endif
126
127! ! Set the inflow values of the dyes, one per segment.
128! ! We know the order: north, south, east, west
129 do m=1,ntr
130 write(name,'("dye_",I2.2)') m
131 write(longname,'("Concentration of dyed_obc Tracer ",I2.2, " on segment ",I2.2)') m, m
132 call tracer_name_lookup(tr_reg, ntr_id, tr_ptr, name)
133
134 do n=1,obc%number_of_segments
135 if (n == m) then
136 dye = 1.0
137 else
138 dye = 0.0
139 endif
140 call register_segment_tracer(tr_ptr, ntr_id, param_file, gv, &
141 obc%segment(n), obc_scalar=dye)
142 enddo
143 enddo
144
146
147!> This subroutine updates the long-channel flow
148subroutine dyed_channel_update_flow(OBC, CS, G, GV, US, h, Time)
149 type(ocean_obc_type), pointer :: obc !< This open boundary condition type specifies
150 !! whether, where, and what open boundary
151 !! conditions are used.
152 type(dyed_channel_obc_cs), pointer :: cs !< Dyed channel control structure.
153 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
154 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
155 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
156 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< layer thickness [H ~> m or kg m-2]
157 type(time_type), intent(in) :: time !< model time.
158
159 ! Local variables
160 real :: flow ! The OBC velocity [L T-1 ~> m s-1]
161 real :: pi ! 3.1415926535... [nondim]
162 real :: time_sec ! The elapsed time since the start of the calendar [T ~> s]
163 real :: fixed_thickness ! A fixed layer thickness, hard-coded to 1 mks unit, that is used to
164 ! reproduce a bug with the older versions of this code [H ~> m or kg m-2]
165 logical :: cross_channel ! True if the segment runs across the channel
166 integer :: turns ! Number of index quarter turns
167 integer :: i, j, k, l_seg, isd, ied, jsd, jed
168 integer :: isdb, iedb, jsdb, jedb, is, ie, js, je
169 type(obc_segment_type), pointer :: segment => null()
170
171 if (.not.associated(obc)) call mom_error(fatal, 'dyed_channel_initialization.F90: '// &
172 'dyed_channel_update_flow() was called but OBC type was not initialized!')
173
174 time_sec = time_to_real(time, scale=us%s_to_T)
175 pi = 4.0*atan(1.0)
176
177 turns = modulo(g%HI%turns, 4)
178
179 do l_seg=1, obc%number_of_segments
180 segment => obc%segment(l_seg)
181 if (.not. segment%on_pe) cycle
182 if (segment%gradient) cycle
183 if (segment%oblique .and. (.not. segment%nudged) .and. (.not. segment%Flather)) cycle
184
185 if (cs%frequency == 0.0) then
186 flow = cs%zonal_flow
187 else
188 flow = cs%zonal_flow + cs%tidal_amp * cos(2 * pi * cs%frequency * time_sec)
189 endif
190 if ((turns==2) .or. (turns==3)) flow = -1.0 * flow
191
192 isd = segment%HI%isd ; ied = segment%HI%ied
193 jsd = segment%HI%jsd ; jed = segment%HI%jed
194 isdb = segment%HI%IsdB ; iedb = segment%HI%IedB
195 jsdb = segment%HI%JsdB ; jedb = segment%HI%JedB
196 if (segment%is_E_or_W) then
197 is = isdb ; ie = iedb ; js = jsd ; je = jed
198 else
199 is = isd ; ie = ied ; js = jsdb ; je = jedb
200 endif
201 cross_channel = ((segment%is_E_or_W .and. ((turns==0) .or. (turns==2))) .or. &
202 (segment%is_N_or_S .and. ((turns==1) .or. (turns==3))))
203
204 if ((segment%specified .or. segment%nudged) .and. cross_channel) then
205 do k=1,gv%ke ; do j=js,je ; do i=is,ie
206 segment%normal_vel(i,j,k) = flow
207 enddo ; enddo ; enddo
208 endif
209
210 if (segment%specified .and. cross_channel) then
211 if (cs%OBC_transport_bug) then
212 fixed_thickness = 1.0 / gv%H_to_mks ! This replicates the prevoius answers without rescaling.
213 if ((segment%direction == obc_direction_w) .or. (segment%direction == obc_direction_e)) then
214 do k=1,gv%ke ; do j=jsd,jed ; do i=isdb,iedb
215 segment%normal_trans(i,j,k) = flow * g%dyCu(i,j) * fixed_thickness
216 enddo ; enddo ; enddo
217 elseif ((segment%direction == obc_direction_s) .or. (segment%direction == obc_direction_n)) then
218 do k=1,gv%ke ; do j=jsdb,jedb ; do i=isd,ied
219 segment%normal_trans(i,j,k) = flow * g%dxCv(i,j) * fixed_thickness
220 enddo ; enddo ; enddo
221 endif
222 else
223 if (segment%direction == obc_direction_w) then
224 do k=1,gv%ke ; do j=jsd,jed ; do i=isdb,iedb
225 segment%normal_trans(i,j,k) = flow * g%dyCu(i,j) * h(i+1,j,k)
226 enddo ; enddo ; enddo
227 elseif (segment%direction == obc_direction_e) then
228 do k=1,gv%ke ; do j=jsd,jed ; do i=isdb,iedb
229 segment%normal_trans(i,j,k) = flow * g%dyCu(i,j) * h(i,j,k)
230 enddo ; enddo ; enddo
231 elseif (segment%direction == obc_direction_s) then
232 do k=1,gv%ke ; do j=jsdb,jedb ; do i=isd,ied
233 segment%normal_trans(i,j,k) = flow * g%dxCv(i,j) * h(i,j+1,k)
234 enddo ; enddo ; enddo
235 elseif (segment%direction == obc_direction_n) then
236 do k=1,gv%ke ; do j=jsdb,jedb ; do i=isd,ied
237 segment%normal_trans(i,j,k) = flow * g%dxCv(i,j) * h(i,j,k)
238 enddo ; enddo ; enddo
239 endif
240 endif
241 endif
242
243 if (cross_channel) then
244 do j=js,je ; do i=is,ie
245 segment%normal_vel_bt(i,j) = flow
246 enddo ; enddo
247 else
248 do j=js,je ; do i=is,ie
249 segment%normal_vel_bt(i,j) = 0.0
250 enddo ; enddo
251 endif
252
253 enddo
254
255end subroutine dyed_channel_update_flow
256
257!> \namespace dyed_channel_initialization
258!!
259!! Setting dyes, one for painting the inflow on each side.