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