MOM_controlled_forcing.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!> Use control-theory to adjust the surface heat flux and precipitation.
6!!
7!! Adjustments are based on the time-mean or periodically (seasonally) varying
8!! anomalies from the observed state.
9!!
10!! The techniques behind this are described in Hallberg and Adcroft (2018, in prep.).
12
14use mom_diag_mediator, only : register_diag_field, diag_ctrl, safe_alloc_ptr
15use mom_domains, only : pass_var, pass_vector, agrid, to_south, to_west, to_all
16use mom_error_handler, only : mom_error, fatal, warning, mom_mesg, is_root_pe
17use mom_file_parser, only : read_param, get_param, log_param, log_version, param_file_type
18use mom_forcing_type, only : forcing
19use mom_grid, only : ocean_grid_type
20use mom_restart, only : register_restart_field, mom_restart_cs
21use mom_time_manager, only : time_type, operator(+), operator(/), operator(-)
22use mom_time_manager, only : get_date, set_date
23use mom_time_manager, only : time_type_to_real, real_to_time
25use mom_variables, only : surface
26
27implicit none ; private
28
29#include <MOM_memory.h>
30
33
34!> Control structure for MOM_controlled_forcing
35type, public :: ctrl_forcing_cs ; private
36 logical :: use_temperature !< If true, temperature and salinity are used as state variables.
37 logical :: do_integrated !< If true, use time-integrated anomalies to control the surface state.
38 integer :: num_cycle !< The number of elements in the forcing cycle.
39 real :: heat_int_rate !< The rate at which heating anomalies accumulate [T-1 ~> s-1]
40 real :: prec_int_rate !< The rate at which precipitation anomalies accumulate [T-1 ~> s-1]
41 real :: heat_cyc_rate !< The rate at which cyclical heating anomalies accumulate [T-1 ~> s-1]
42 real :: prec_cyc_rate !< The rate at which cyclical precipitation anomalies
43 !! accumulate [T-1 ~> s-1]
44 real :: len2 !< The square of the length scale over which the anomalies
45 !! are smoothed via a Laplacian filter [L2 ~> m2]
46 real :: lam_heat !< A constant of proportionality between SST anomalies
47 !! and heat fluxes [Q R Z T-1 C-1 ~> W m-2 degC-1]
48 real :: lam_prec !< A constant of proportionality between SSS anomalies
49 !! (normalised by mean SSS) and precipitation [R Z T-1 ~> kg m-2 s-1]
50 real :: lam_cyc_heat !< A constant of proportionality between cyclical SST
51 !! anomalies and corrective heat fluxes [Q R Z T-1 C-1 ~> W m-2 degC-1]
52 real :: lam_cyc_prec !< A constant of proportionality between cyclical SSS
53 !! anomalies (normalised by mean SSS) and corrective
54 !! precipitation [R Z T-1 ~> kg m-2 s-1]
55
56 real, pointer, dimension(:,:) :: &
57 heat_0 => null(), & !< The non-periodic integrative corrective heat flux that has been
58 !! evolved to control mean SST anomalies [Q R Z T-1 ~> W m-2]
59 precip_0 => null() !< The non-periodic integrative corrective precipitation that has been
60 !! evolved to control mean SSS anomalies [R Z T-1 ~> kg m-2 s-1]
61
62 ! The final dimension of each of the six variables that follow is for the periodic bins.
63 real, pointer, dimension(:,:,:) :: &
64 heat_cyc => null(), & !< The periodic integrative corrective heat flux that has been evolved
65 !! to control periodic (seasonal) SST anomalies [Q R Z T-1 ~> W m-2].
66 !! The third dimension is the periodic bins.
67 precip_cyc => null() !< The non-periodic integrative corrective precipitation that has been
68 !! evolved to control periodic (seasonal) SSS anomalies [R Z T-1 ~> kg m-2 s-1].
69 !! The third dimension is the periodic bins.
70 real, pointer, dimension(:) :: &
71 avg_time => null() !< The accumulated averaging time in each part of the cycle [T ~> s] or
72 !! a negative value to indicate that the variables like avg_SST_anom are
73 !! the actual averages, and not time integrals.
74 !! The dimension is the periodic bins.
75 real, pointer, dimension(:,:,:) :: &
76 avg_sst_anom => null(), & !< The time-averaged periodic sea surface temperature anomalies [C ~> degC],
77 !! or (at some points in the code), the time-integrated periodic
78 !! temperature anomalies [T C ~> s degC].
79 !! The third dimension is the periodic bins.
80 avg_sss_anom => null(), & !< The time-averaged periodic sea surface salinity anomalies [S ~> ppt],
81 !! or (at some points in the code), the time-integrated periodic
82 !! salinity anomalies [T S ~> s ppt].
83 !! The third dimension is the periodic bins.
84 avg_sss => null() !< The time-averaged periodic sea surface salinities [S ~> ppt], or (at
85 !! some points in the code), the time-integrated periodic
86 !! salinities [T S ~> s ppt].
87 !! The third dimension is the periodic bins.
88
89 type(diag_ctrl), pointer :: diag => null() !< A structure that is used to
90 !! regulate the timing of diagnostic output.
91 integer :: id_heat_0 = -1 !< Diagnostic handle for the steady heat flux
92 integer :: id_prec_0 = -1 !< Diagnostic handle for the steady precipitation
93end type ctrl_forcing_cs
94
95contains
96
97!> This subroutine determines corrective surface forcing fields using simple control theory.
98subroutine apply_ctrl_forcing(SST_anom, SSS_anom, SSS_mean, virt_heat, virt_precip, &
99 day_start, dt, G, US, CS)
100 type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure
101 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: sst_anom !< The sea surface temperature anomalies [C ~> degC]
102 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: sss_anom !< The sea surface salinity anomlies [S ~> ppt]
103 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: sss_mean !< The mean sea surface salinity [S ~> ppt]
104 real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: virt_heat !< Virtual (corrective) heat
105 !! fluxes that are augmented in this
106 !! subroutine [Q R Z T-1 ~> W m-2]
107 real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: virt_precip !< Virtual (corrective)
108 !! precipitation fluxes that are augmented
109 !! in this subroutine [R Z T-1 ~> kg m-2 s-1]
110 type(time_type), intent(in) :: day_start !< Start time of the fluxes.
111 real, intent(in) :: dt !< Length of time over which these fluxes
112 !! will be applied [T ~> s]
113 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
114 type(ctrl_forcing_cs), pointer :: cs !< A pointer to the control structure returned
115 !! by a previous call to ctrl_forcing_init.
116
117 ! Local variables
118 real, dimension(SZIB_(G),SZJ_(G)) :: &
119 flux_heat_x, & ! Zonal smoothing flux of the virtual heat fluxes [L2 Q R Z T-1 ~> W]
120 flux_prec_x ! Zonal smoothing flux of the virtual precipitation [L2 R Z T-1 ~> kg s-1]
121 real, dimension(SZI_(G),SZJB_(G)) :: &
122 flux_heat_y, & ! Meridional smoothing flux of the virtual heat fluxes [L2 Q R Z T-1 ~> W]
123 flux_prec_y ! Meridional smoothing flux of the virtual precipitation [L2 R Z T-1 ~> kg s-1]
124 type(time_type) :: day_end
125 real :: coef ! A heat-flux coefficient [L2 ~> m2]
126 real :: mr_st, mr_end, mr_mid ! Position of various times in the periodic cycle [nondim]
127 real :: mr_prev, mr_next ! Position of various times in the periodic cycle [nondim]
128 real :: dt_wt ! The timestep times a fractional weight used to accumulate averages [T ~> s]
129 real :: dt_heat_rate, dt_prec_rate ! Timestep times the flux accumulation rate [nondim]
130 real :: dt1_heat_rate, dt1_prec_rate, dt2_heat_rate, dt2_prec_rate ! [nondim]
131 real :: wt_per1, wt_st, wt_end, wt_mid ! Averaging weights [nondim]
132 integer :: m_st, m_end, m_mid, m_u1, m_u2, m_u3 ! Indices (nominally months) in the periodic cycle
133 integer :: yr, mon, day, hr, min, sec
134 integer :: i, j, is, ie, js, je
135
136 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
137
138 if (.not.associated(cs)) return
139 if ((cs%num_cycle <= 0) .and. (.not.cs%do_integrated)) return
140
141 day_end = day_start + real_to_time(us%T_to_s*dt)
142
143 do j=js,je ; do i=is,ie
144 virt_heat(i,j) = 0.0 ; virt_precip(i,j) = 0.0
145 enddo ; enddo
146
147 if (cs%do_integrated) then
148 dt_heat_rate = dt * cs%heat_int_rate
149 dt_prec_rate = dt * cs%prec_int_rate
150 call pass_var(cs%heat_0, g%Domain, complete=.false.)
151 call pass_var(cs%precip_0, g%Domain)
152
153 do j=js,je ; do i=is-1,ie
154 coef = cs%Len2 * (g%dy_Cu(i,j)*g%IdxCu(i,j))
155 flux_heat_x(i,j) = coef * (cs%heat_0(i,j) - cs%heat_0(i+1,j))
156 flux_prec_x(i,j) = coef * (cs%precip_0(i,j) - cs%precip_0(i+1,j))
157 enddo ; enddo
158 do j=js-1,je ; do i=is,ie
159 coef = cs%Len2 * (g%dx_Cv(i,j)*g%IdyCv(i,j))
160 flux_heat_y(i,j) = coef * (cs%heat_0(i,j) - cs%heat_0(i,j+1))
161 flux_prec_y(i,j) = coef * (cs%precip_0(i,j) - cs%precip_0(i,j+1))
162 enddo ; enddo
163 do j=js,je ; do i=is,ie
164 cs%heat_0(i,j) = cs%heat_0(i,j) + dt_heat_rate * ( &
165 -cs%lam_heat*g%mask2dT(i,j)*sst_anom(i,j) + &
166 (g%IareaT(i,j) * ((flux_heat_x(i-1,j) - flux_heat_x(i,j)) + &
167 (flux_heat_y(i,j-1) - flux_heat_y(i,j))) ) )
168
169 cs%precip_0(i,j) = cs%precip_0(i,j) + dt_prec_rate * ( &
170 cs%lam_prec * g%mask2dT(i,j)*(sss_anom(i,j) / sss_mean(i,j)) + &
171 (g%IareaT(i,j) * ((flux_prec_x(i-1,j) - flux_prec_x(i,j)) + &
172 (flux_prec_y(i,j-1) - flux_prec_y(i,j))) ) )
173
174 virt_heat(i,j) = virt_heat(i,j) + cs%heat_0(i,j)
175 virt_precip(i,j) = virt_precip(i,j) + cs%precip_0(i,j)
176 enddo ; enddo
177 endif
178
179 if (cs%num_cycle > 0) then
180 ! Determine the current period, with values that run from 0 to CS%num_cycle.
181 call get_date(day_start, yr, mon, day, hr, min, sec)
182 mr_st = cs%num_cycle * (time_type_to_real(day_start - set_date(yr, 1, 1)) / &
183 time_type_to_real(set_date(yr+1, 1, 1) - set_date(yr, 1, 1)))
184
185 call get_date(day_end, yr, mon, day, hr, min, sec)
186 mr_end = cs%num_cycle * (time_type_to_real(day_end - set_date(yr, 1, 1)) / &
187 time_type_to_real(set_date(yr+1, 1, 1) - set_date(yr, 1, 1)))
188
189 ! The Chapeau functions are centered at whole integer values that are nominally
190 ! the end of the month to enable simple conversion from the fractional-years times
191 ! CS%num_cycle.
192
193 ! The month-average temperatures have as an index the month number.
194
195 m_end = periodic_int(real(ceiling(mr_end)), cs%num_cycle)
196 m_mid = periodic_int(real(ceiling(mr_st)), cs%num_cycle)
197 m_st = periodic_int(mr_st, cs%num_cycle)
198
199 mr_st = periodic_real(mr_st, cs%num_cycle)
200 mr_end = periodic_real(mr_end, cs%num_cycle)
201 ! mr_mid = periodic_real(ceiling(mr_st), CS%num_cycle)
202 mr_prev = periodic_real(real(floor(mr_st)), cs%num_cycle)
203 mr_next = periodic_real(real(m_end), cs%num_cycle)
204 if (m_mid == m_end) then ; mr_mid = mr_end ! There is only one cell.
205 else ; mr_mid = periodic_real(real(m_mid), cs%num_cycle) ; endif
206
207 ! There may be two cells that run from mr_st to mr_mid and mr_mid to mr_end.
208
209 ! The values of m for weights are all calculated relative to mr_prev, so
210 ! check whether mr_mid, etc., need to be shifted by CS%num_cycle, so that these
211 ! values satisfiy mr_prev <= mr_st < mr_mid <= mr_end <= mr_next.
212 if (mr_st < mr_prev) mr_prev = mr_prev - cs%num_cycle
213 if (mr_mid < mr_st) mr_mid = mr_mid + cs%num_cycle
214 if (mr_end < mr_st) mr_end = mr_end + cs%num_cycle
215 if (mr_next < mr_prev) mr_next = mr_next + cs%num_cycle
216
217 !### These might be removed later - they are to check the coding.
218 if ((mr_mid < mr_st) .or. (mr_mid > mr_prev + 1.)) call mom_error(fatal, &
219 "apply ctrl_forcing: m_mid interpolation out of bounds; fix the code.")
220 if ((mr_end < mr_st) .or. (mr_end > mr_prev + 2.)) call mom_error(fatal, &
221 "apply ctrl_forcing: m_end interpolation out of bounds; fix the code.")
222 if (mr_end > mr_next) call mom_error(fatal, &
223 "apply ctrl_forcing: mr_next interpolation out of bounds; fix the code.")
224
225 wt_per1 = 1.0
226 if (mr_mid < mr_end) wt_per1 = (mr_mid - mr_st) / (mr_end - mr_st)
227
228 ! Find the 3 Chapeau-function weights, bearing in mind that m_end may be m_mid.
229 wt_st = wt_per1 * (1. + (mr_prev - 0.5*(mr_st + mr_mid)))
230 wt_end = (1.0-wt_per1) * (1. + (0.5*(mr_end + mr_mid) - mr_next))
231 wt_mid = 1.0 - (wt_st + wt_end)
232 if ((wt_st < 0.0) .or. (wt_end < 0.0) .or. (wt_mid < 0.0)) &
233 call mom_error(fatal, "apply_ctrl_forcing: Negative m weights")
234 if ((wt_st > 1.0) .or. (wt_end > 1.0) .or. (wt_mid > 1.0)) &
235 call mom_error(fatal, "apply_ctrl_forcing: Excessive m weights")
236
237 ! Add to vert_heat and vert_precip.
238 do j=js,je ; do i=is,ie
239 virt_heat(i,j) = virt_heat(i,j) + (wt_st * cs%heat_cyc(i,j,m_st) + &
240 (wt_mid * cs%heat_cyc(i,j,m_mid) + &
241 wt_end * cs%heat_cyc(i,j,m_end)))
242 virt_precip(i,j) = virt_precip(i,j) + (wt_st * cs%precip_cyc(i,j,m_st) + &
243 (wt_mid * cs%precip_cyc(i,j,m_mid) + &
244 wt_end * cs%precip_cyc(i,j,m_end)))
245 enddo ; enddo
246
247 ! If different from the last period, take the average and determine the
248 ! chapeau weighting
249
250 ! The Chapeau functions are centered at whole integer values that are nominally
251 ! the end of the month to enable simple conversion from the fractional-years times
252 ! CS%num_cycle.
253
254 ! The month-average temperatures have as an index the month number, so the averages
255 ! apply to indicies m_end and m_mid.
256
257 if (cs%avg_time(m_end) <= 0.0) then ! zero out the averages.
258 cs%avg_time(m_end) = 0.0
259 do j=js,je ; do i=is,ie
260 cs%avg_SST_anom(i,j,m_end) = 0.0
261 cs%avg_SSS_anom(i,j,m_end) = 0.0 ; cs%avg_SSS(i,j,m_end) = 0.0
262 enddo ; enddo
263 endif
264 if (cs%avg_time(m_mid) <= 0.0) then ! zero out the averages.
265 cs%avg_time(m_mid) = 0.0
266 do j=js,je ; do i=is,ie
267 cs%avg_SST_anom(i,j,m_mid) = 0.0
268 cs%avg_SSS_anom(i,j,m_mid) = 0.0 ; cs%avg_SSS(i,j,m_mid) = 0.0
269 enddo ; enddo
270 endif
271
272 ! Accumulate the average anomalies for this period.
273 dt_wt = wt_per1 * dt
274 cs%avg_time(m_mid) = cs%avg_time(m_mid) + dt_wt
275 ! These loops temporarily change the units of the CS%avg_ variables to [C T ~> degC s]
276 ! or [S T ~> ppt s].
277 do j=js,je ; do i=is,ie
278 cs%avg_SST_anom(i,j,m_mid) = cs%avg_SST_anom(i,j,m_mid) + &
279 dt_wt * g%mask2dT(i,j) * sst_anom(i,j)
280 cs%avg_SSS_anom(i,j,m_mid) = cs%avg_SSS_anom(i,j,m_mid) + &
281 dt_wt * g%mask2dT(i,j) * sss_anom(i,j)
282 cs%avg_SSS(i,j,m_mid) = cs%avg_SSS(i,j,m_mid) + dt_wt * sss_mean(i,j)
283 enddo ; enddo
284 if (wt_per1 < 1.0) then
285 dt_wt = (1.0-wt_per1) * dt
286 cs%avg_time(m_end) = cs%avg_time(m_end) + dt_wt
287 do j=js,je ; do i=is,ie
288 cs%avg_SST_anom(i,j,m_end) = cs%avg_SST_anom(i,j,m_end) + &
289 dt_wt * g%mask2dT(i,j) * sst_anom(i,j)
290 cs%avg_SSS_anom(i,j,m_end) = cs%avg_SSS_anom(i,j,m_end) + &
291 dt_wt * g%mask2dT(i,j) * sss_anom(i,j)
292 cs%avg_SSS(i,j,m_end) = cs%avg_SSS(i,j,m_end) + dt_wt * sss_mean(i,j)
293 enddo ; enddo
294 endif
295
296 ! Update the Chapeau magnitudes for 4 cycles ago.
297 m_u1 = periodic_int(m_st - 4.0, cs%num_cycle)
298 m_u2 = periodic_int(m_st - 3.0, cs%num_cycle)
299 m_u3 = periodic_int(m_st - 2.0, cs%num_cycle)
300
301 ! These loops restore the units of the CS%avg variables to [C ~> degC] or [S ~> ppt]
302 if (cs%avg_time(m_u1) > 0.0) then
303 do j=js,je ; do i=is,ie
304 cs%avg_SST_anom(i,j,m_u1) = cs%avg_SST_anom(i,j,m_u1) / cs%avg_time(m_u1)
305 cs%avg_SSS_anom(i,j,m_u1) = cs%avg_SSS_anom(i,j,m_u1) / cs%avg_time(m_u1)
306 cs%avg_SSS(i,j,m_u1) = cs%avg_SSS(i,j,m_u1) / cs%avg_time(m_u1)
307 enddo ; enddo
308 cs%avg_time(m_u1) = -1.0
309 endif
310 if (cs%avg_time(m_u2) > 0.0) then
311 do j=js,je ; do i=is,ie
312 cs%avg_SST_anom(i,j,m_u2) = cs%avg_SST_anom(i,j,m_u2) / cs%avg_time(m_u2)
313 cs%avg_SSS_anom(i,j,m_u2) = cs%avg_SSS_anom(i,j,m_u2) / cs%avg_time(m_u2)
314 cs%avg_SSS(i,j,m_u2) = cs%avg_SSS(i,j,m_u2) / cs%avg_time(m_u2)
315 enddo ; enddo
316 cs%avg_time(m_u2) = -1.0
317 endif
318 if (cs%avg_time(m_u3) > 0.0) then
319 do j=js,je ; do i=is,ie
320 cs%avg_SST_anom(i,j,m_u3) = cs%avg_SST_anom(i,j,m_u3) / cs%avg_time(m_u3)
321 cs%avg_SSS_anom(i,j,m_u3) = cs%avg_SSS_anom(i,j,m_u3) / cs%avg_time(m_u3)
322 cs%avg_SSS(i,j,m_u3) = cs%avg_SSS(i,j,m_u3) / cs%avg_time(m_u3)
323 enddo ; enddo
324 cs%avg_time(m_u3) = -1.0
325 endif
326
327 dt1_heat_rate = wt_per1 * dt * cs%heat_cyc_rate
328 dt1_prec_rate = wt_per1 * dt * cs%prec_cyc_rate
329 dt2_heat_rate = (1.0-wt_per1) * dt * cs%heat_cyc_rate
330 dt2_prec_rate = (1.0-wt_per1) * dt * cs%prec_cyc_rate
331
332 if (wt_per1 < 1.0) then
333 call pass_var(cs%heat_cyc(:,:,m_u2), g%Domain, complete=.false.)
334 call pass_var(cs%precip_cyc(:,:,m_u2), g%Domain, complete=.false.)
335 endif
336 call pass_var(cs%heat_cyc(:,:,m_u1), g%Domain, complete=.false.)
337 call pass_var(cs%precip_cyc(:,:,m_u1), g%Domain)
338
339 if ((cs%avg_time(m_u1) == -1.0) .and. (cs%avg_time(m_u2) == -1.0)) then
340 do j=js,je ; do i=is-1,ie
341 coef = cs%Len2 * (g%dy_Cu(i,j)*g%IdxCu(i,j))
342 flux_heat_x(i,j) = coef * (cs%heat_cyc(i,j,m_u1) - cs%heat_cyc(i+1,j,m_u1))
343 flux_prec_x(i,j) = coef * (cs%precip_cyc(i,j,m_u1) - cs%precip_cyc(i+1,j,m_u1))
344 enddo ; enddo
345 do j=js-1,je ; do i=is,ie
346 coef = cs%Len2 * (g%dx_Cv(i,j)*g%IdyCv(i,j))
347 flux_heat_y(i,j) = coef * (cs%heat_cyc(i,j,m_u1) - cs%heat_cyc(i,j+1,m_u1))
348 flux_prec_y(i,j) = coef * (cs%precip_cyc(i,j,m_u1) - cs%precip_cyc(i,j+1,m_u1))
349 enddo ; enddo
350 do j=js,je ; do i=is,ie
351 cs%heat_cyc(i,j,m_u1) = cs%heat_cyc(i,j,m_u1) + dt1_heat_rate * ( &
352 -cs%lam_cyc_heat*(cs%avg_SST_anom(i,j,m_u2) - cs%avg_SST_anom(i,j,m_u1)) + &
353 (g%IareaT(i,j) * ((flux_heat_x(i-1,j) - flux_heat_x(i,j)) + &
354 (flux_heat_y(i,j-1) - flux_heat_y(i,j))) ) )
355
356 cs%precip_cyc(i,j,m_u1) = cs%precip_cyc(i,j,m_u1) + dt1_prec_rate * ( &
357 cs%lam_prec * (cs%avg_SSS_anom(i,j,m_u2) - cs%avg_SSS_anom(i,j,m_u1)) / &
358 (0.5*(cs%avg_SSS(i,j,m_u2) + cs%avg_SSS(i,j,m_u1))) + &
359 (g%IareaT(i,j) * ((flux_prec_x(i-1,j) - flux_prec_x(i,j)) + &
360 (flux_prec_y(i,j-1) - flux_prec_y(i,j))) ) )
361 enddo ; enddo
362 endif
363
364 if ((wt_per1 < 1.0) .and. (cs%avg_time(m_u1) == -1.0) .and. (cs%avg_time(m_u2) == -1.0)) then
365 do j=js,je ; do i=is-1,ie
366 coef = cs%Len2 * (g%dy_Cu(i,j)*g%IdxCu(i,j))
367 flux_heat_x(i,j) = coef * (cs%heat_cyc(i,j,m_u2) - cs%heat_cyc(i+1,j,m_u2))
368 flux_prec_x(i,j) = coef * (cs%precip_cyc(i,j,m_u2) - cs%precip_cyc(i+1,j,m_u2))
369 enddo ; enddo
370 do j=js-1,je ; do i=is,ie
371 coef = cs%Len2 * (g%dx_Cv(i,j)*g%IdyCv(i,j))
372 flux_heat_y(i,j) = coef * (cs%heat_cyc(i,j,m_u2) - cs%heat_cyc(i,j+1,m_u2))
373 flux_prec_y(i,j) = coef * (cs%precip_cyc(i,j,m_u2) - cs%precip_cyc(i,j+1,m_u2))
374 enddo ; enddo
375 do j=js,je ; do i=is,ie
376 cs%heat_cyc(i,j,m_u2) = cs%heat_cyc(i,j,m_u2) + dt1_heat_rate * ( &
377 -cs%lam_cyc_heat*(cs%avg_SST_anom(i,j,m_u3) - cs%avg_SST_anom(i,j,m_u2)) + &
378 (g%IareaT(i,j) * ((flux_heat_x(i-1,j) - flux_heat_x(i,j)) + &
379 (flux_heat_y(i,j-1) - flux_heat_y(i,j))) ) )
380
381 cs%precip_cyc(i,j,m_u2) = cs%precip_cyc(i,j,m_u2) + dt1_prec_rate * ( &
382 cs%lam_prec * (cs%avg_SSS_anom(i,j,m_u3) - cs%avg_SSS_anom(i,j,m_u2)) / &
383 (0.5*(cs%avg_SSS(i,j,m_u3) + cs%avg_SSS(i,j,m_u2))) + &
384 (g%IareaT(i,j) * ((flux_prec_x(i-1,j) - flux_prec_x(i,j)) + &
385 (flux_prec_y(i,j-1) - flux_prec_y(i,j))) ) )
386 enddo ; enddo
387 endif
388
389 endif ! (CS%num_cycle > 0)
390
391 if (cs%do_integrated .and. ((cs%id_heat_0 > 0) .or. (cs%id_prec_0 > 0))) then
392 call enable_averages(dt, day_start + real_to_time(us%T_to_s*dt), cs%diag)
393 if (cs%id_heat_0 > 0) call post_data(cs%id_heat_0, cs%heat_0, cs%diag)
394 if (cs%id_prec_0 > 0) call post_data(cs%id_prec_0, cs%precip_0, cs%diag)
395 call disable_averaging(cs%diag)
396 endif
397
398end subroutine apply_ctrl_forcing
399
400!> This function maps rval into an integer in the range from 1 to num_period.
401function periodic_int(rval, num_period) result (m)
402 real, intent(in) :: rval !< Input for mapping [nondim]
403 integer, intent(in) :: num_period !< Maximum output.
404 integer :: m !< Return value.
405
406 m = floor(rval)
407 if (m <= 0) then
408 m = m + num_period * (1 + (abs(m) / num_period))
409 elseif (m > num_period) then
410 m = m - num_period * ((m-1) / num_period)
411 endif
412end function
413
414!> This function shifts rval by an integer multiple of num_period so that
415!! 0 <= val_out < num_period.
416function periodic_real(rval, num_period) result(val_out)
417 real, intent(in) :: rval !< Input to be shifted into valid range [nondim]
418 integer, intent(in) :: num_period !< Maximum valid value.
419 real :: val_out !< Return value [nondim]
420 integer :: nshft
421
422 if (rval < 0) then ; nshft = floor(abs(rval) / num_period) + 1
423 elseif (rval < num_period) then ; nshft = 0
424 else ; nshft = -1*floor(rval / num_period) ; endif
425
426 val_out = rval + nshft * num_period
427end function
428
429
430!> This subroutine is used to allocate and register any fields in this module
431!! that should be written to or read from the restart file.
432subroutine register_ctrl_forcing_restarts(G, US, param_file, CS, restart_CS)
433 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
434 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
435 type(param_file_type), intent(in) :: param_file !< A structure indicating the
436 !! open file to parse for model
437 !! parameter values.
438 type(ctrl_forcing_cs), pointer :: cs !< A pointer that is set to point to the
439 !! control structure for this module.
440 type(mom_restart_cs), intent(inout) :: restart_cs !< MOM restart control struct
441
442 logical :: controlled, use_temperature
443 character (len=8) :: period_str
444 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
445 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
446 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
447
448 if (associated(cs)) then
449 call mom_error(warning, "register_ctrl_forcing_restarts called "//&
450 "with an associated control structure.")
451 return
452 endif
453
454 controlled = .false.
455 call read_param(param_file, "CONTROLLED_FORCING", controlled)
456 if (.not.controlled) return
457
458 use_temperature = .true.
459 call read_param(param_file, "ENABLE_THERMODYNAMICS", use_temperature)
460 if (.not.use_temperature) call mom_error(fatal, &
461 "register_ctrl_forcing_restarts: CONTROLLED_FORCING only works with "//&
462 "ENABLE_THERMODYNAMICS defined.")
463
464 allocate(cs)
465
466 cs%do_integrated = .true. ; cs%num_cycle = 0
467 call read_param(param_file, "CTRL_FORCE_INTEGRATED", cs%do_integrated)
468 call read_param(param_file, "CTRL_FORCE_NUM_CYCLE", cs%num_cycle)
469
470 if (cs%do_integrated) then
471 allocate(cs%heat_0(isd:ied,jsd:jed), source=0.0)
472 allocate(cs%precip_0(isd:ied,jsd:jed), source=0.0)
473
474 call register_restart_field(cs%heat_0, "Ctrl_heat", .false., restart_cs, &
475 longname="Control Integrative Heating", &
476 units="W m-2", conversion=us%QRZ_T_to_W_m2, z_grid='1')
477 call register_restart_field(cs%precip_0, "Ctrl_precip", .false., restart_cs, &
478 longname="Control Integrative Precipitation", &
479 units="kg m-2 s-1", conversion=us%RZ_T_to_kg_m2s, z_grid='1')
480 endif
481
482 if (cs%num_cycle > 0) then
483 allocate(cs%heat_cyc(isd:ied,jsd:jed,cs%num_cycle), source=0.0)
484 allocate(cs%precip_cyc(isd:ied,jsd:jed,cs%num_cycle), source=0.0)
485 allocate(cs%avg_time(cs%num_cycle), source=0.0)
486 allocate(cs%avg_SST_anom(isd:ied,jsd:jed,cs%num_cycle), source=0.0)
487 allocate(cs%avg_SSS_anom(isd:ied,jsd:jed,cs%num_cycle), source=0.0)
488 allocate(cs%avg_SSS(isd:ied,jsd:jed,cs%num_cycle), source=0.0)
489
490 write (period_str, '("p ",I0)') cs%num_cycle
491
492 call register_restart_field(cs%heat_cyc, "Ctrl_heat_cycle", .false., restart_cs, &
493 longname="Cyclical Control Heating", &
494 units="W m-2", conversion=us%QRZ_T_to_W_m2, z_grid='1', t_grid=period_str)
495 call register_restart_field(cs%precip_cyc, "Ctrl_precip_cycle", .false., restart_cs, &
496 longname="Cyclical Control Precipitation", &
497 units="kg m-2 s-1", conversion=us%RZ_T_to_kg_m2s, z_grid='1', t_grid=period_str)
498 call register_restart_field(cs%avg_time, "avg_time", .false., restart_cs, &
499 longname="Cyclical accumulated averaging time", &
500 units="sec", conversion=us%T_to_s, z_grid='1', t_grid=period_str)
501 call register_restart_field(cs%avg_SST_anom, "avg_SST_anom", .false., restart_cs, &
502 longname="Cyclical average SST Anomaly", &
503 units="degC", conversion=us%C_to_degC, z_grid='1', t_grid=period_str)
504 call register_restart_field(cs%avg_SSS_anom, "avg_SSS_anom", .false., restart_cs, &
505 longname="Cyclical average SSS Anomaly", &
506 units="g kg-1", conversion=us%S_to_ppt, z_grid='1', t_grid=period_str)
507 call register_restart_field(cs%avg_SSS_anom, "avg_SSS", .false., restart_cs, &
508 longname="Cyclical average SSS", &
509 units="g kg-1", conversion=us%S_to_ppt, z_grid='1', t_grid=period_str)
510 endif
511
513
514!> Set up this modules control structure.
515subroutine controlled_forcing_init(Time, G, US, param_file, diag, CS)
516 type(time_type), intent(in) :: time !< The current model time.
517 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
518 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
519 type(param_file_type), intent(in) :: param_file !< A structure indicating the
520 !! open file to parse for model
521 !! parameter values.
522 type(diag_ctrl), target, intent(in) :: diag !< A structure that is used to regulate
523 !! diagnostic output.
524 type(ctrl_forcing_cs), pointer :: cs !< A pointer that is set to point to the
525 !! control structure for this module.
526
527 ! Local variables
528 real :: smooth_len ! A smoothing lengthscale [L ~> m]
529 logical :: do_integrated
530 integer :: num_cycle
531 integer :: i, j, isc, iec, jsc, jec, m
532 ! This include declares and sets the variable "version".
533# include "version_variable.h"
534 character(len=40) :: mdl = "MOM_controlled_forcing" ! This module's name.
535
536 isc = g%isc ; iec = g%iec ; jsc = g%jsc ; jec = g%jec
537
538 ! These should have already been called.
539 ! call read_param(param_file, "CTRL_FORCE_INTEGRATED", CS%do_integrated)
540 ! call read_param(param_file, "CTRL_FORCE_NUM_CYCLE", CS%num_cycle)
541
542 if (associated(cs)) then
543 do_integrated = cs%do_integrated ; num_cycle = cs%num_cycle
544 else
545 do_integrated = .false. ; num_cycle = 0
546 endif
547
548 ! Read all relevant parameters and write them to the model log.
549 call log_version(param_file, mdl, version, "")
550 call log_param(param_file, mdl, "CTRL_FORCE_INTEGRATED", do_integrated, &
551 "If true, use a PI controller to determine the surface "//&
552 "forcing that is consistent with the observed mean properties.", &
553 default=.false.)
554 call log_param(param_file, mdl, "CTRL_FORCE_NUM_CYCLE", num_cycle, &
555 "The number of cycles per year in the controlled forcing, "//&
556 "or 0 for no cyclic forcing.", default=0)
557
558 if (.not.associated(cs)) return
559
560 cs%diag => diag
561
562 call get_param(param_file, mdl, "CTRL_FORCE_HEAT_INT_RATE", cs%heat_int_rate, &
563 "The integrated rate at which heat flux anomalies are accumulated.", &
564 units="s-1", default=0.0, scale=us%T_to_s)
565 call get_param(param_file, mdl, "CTRL_FORCE_PREC_INT_RATE", cs%prec_int_rate, &
566 "The integrated rate at which precipitation anomalies are accumulated.", &
567 units="s-1", default=0.0, scale=us%T_to_s)
568 call get_param(param_file, mdl, "CTRL_FORCE_HEAT_CYC_RATE", cs%heat_cyc_rate, &
569 "The integrated rate at which cyclical heat flux anomalies are accumulated.", &
570 units="s-1", default=0.0, scale=us%T_to_s)
571 call get_param(param_file, mdl, "CTRL_FORCE_PREC_CYC_RATE", cs%prec_cyc_rate, &
572 "The integrated rate at which cyclical precipitation anomalies are accumulated.", &
573 units="s-1", default=0.0, scale=us%T_to_s)
574 call get_param(param_file, mdl, "CTRL_FORCE_SMOOTH_LENGTH", smooth_len, &
575 "The length scales over which controlled forcing anomalies are smoothed.", &
576 units="m", default=0.0, scale=us%m_to_L)
577 call get_param(param_file, mdl, "CTRL_FORCE_LAMDA_HEAT", cs%lam_heat, &
578 "A constant of proportionality between SST anomalies "//&
579 "and controlling heat fluxes", &
580 units="W m-2 K-1", default=0.0, scale=us%W_m2_to_QRZ_T*us%C_to_degC)
581 call get_param(param_file, mdl, "CTRL_FORCE_LAMDA_PREC", cs%lam_prec, &
582 "A constant of proportionality between SSS anomalies "//&
583 "(normalised by mean SSS) and controlling precipitation.", &
584 units="kg m-2 s-1", default=0.0, scale=us%kg_m2s_to_RZ_T)
585 call get_param(param_file, mdl, "CTRL_FORCE_LAMDA_CYC_HEAT", cs%lam_cyc_heat, &
586 "A constant of proportionality between SST anomalies "//&
587 "and cyclical controlling heat fluxes", &
588 units="W m-2 K-1", default=0.0, scale=us%W_m2_to_QRZ_T*us%C_to_degC)
589 call get_param(param_file, mdl, "CTRL_FORCE_LAMDA_CYC_PREC", cs%lam_cyc_prec, &
590 "A constant of proportionality between SSS anomalies "//&
591 "(normalised by mean SSS) and cyclical controlling precipitation.", &
592 units="kg m-2 s-1", default=0.0, scale=us%kg_m2s_to_RZ_T)
593
594 cs%Len2 = smooth_len**2
595
596 if (cs%do_integrated) then
597 cs%id_heat_0 = register_diag_field('ocean_model', 'Ctrl_heat', diag%axesT1, time, &
598 'Control Corrective Heating', 'W m-2', conversion=us%QRZ_T_to_W_m2)
599 cs%id_prec_0 = register_diag_field('ocean_model', 'Ctrl_prec', diag%axesT1, time, &
600 'Control Corrective Precipitation', 'kg m-2 s-1', conversion=us%RZ_T_to_kg_m2s)
601 endif
602
603end subroutine controlled_forcing_init
604
605!> Clean up this modules control structure.
606subroutine controlled_forcing_end(CS)
607 type(ctrl_forcing_cs), pointer :: cs !< A pointer to the control structure
608 !! returned by a previous call to
609 !! controlled_forcing_init, it will be
610 !! deallocated here.
611
612 if (associated(cs)) then
613 if (associated(cs%heat_0)) deallocate(cs%heat_0)
614 if (associated(cs%precip_0)) deallocate(cs%precip_0)
615 if (associated(cs%heat_cyc)) deallocate(cs%heat_cyc)
616 if (associated(cs%precip_cyc)) deallocate(cs%precip_cyc)
617 if (associated(cs%avg_SST_anom)) deallocate(cs%avg_SST_anom)
618 if (associated(cs%avg_SSS_anom)) deallocate(cs%avg_SSS_anom)
619 if (associated(cs%avg_SSS)) deallocate(cs%avg_SSS)
620
621 deallocate(cs)
622 endif
623 cs => null()
624
625end subroutine controlled_forcing_end
626
627!> \namespace mom_controlled_forcing
628!! *
629!! By Robert Hallberg, July 2011 *
630!! *
631!! This program contains the subroutines that use control-theory *
632!! to adjust the surface heat flux and precipitation, based on the *
633!! time-mean or periodically (seasonally) varying anomalies from the *
634!! observed state. The techniques behind this are described in *
635!! Hallberg and Adcroft (2011, in prep.). *
636!! *
637!! Macros written all in capital letters are defined in MOM_memory.h. *
638!! *
639!! A small fragment of the grid is shown below: *
640!! *
641!! j+1 x ^ x ^ x At x: q *
642!! j+1 > o > o > At ^: v, tauy *
643!! j x ^ x ^ x At >: u, taux *
644!! j > o > o > At o: h, fluxes. *
645!! j-1 x ^ x ^ x *
646!! i-1 i i+1 At x & ^: *
647!! i i+1 At > & o: *
648!! *
649!! The boundaries always run through q grid points (x). *
650end module mom_controlled_forcing