MOM_oda_incupd.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!> This module contains the routines used to apply incremental updates
6!! from data assimilation.
7!
8!! Applying incremental updates requires the following:
9!! 1. initialize_oda_incupd_fixed and initialize_oda_incupd
10!! 2. set_up_oda_incupd_field (tracers) and set_up_oda_incupd_vel_field (vel)
11!! 3. calc_oda_increments (if using full fields input)
12!! 4. apply_oda_incupd
13!! 5. output_oda_incupd_inc (output increment if using full fields input)
14!! 6. init_oda_incupd_diags (to output increments in diagnostics)
15!! 7. oda_incupd_end (not being used for now)
16
17module mom_oda_incupd
18
19
20use mom_array_transform, only : rotate_array
21use mom_coms, only : sum_across_pes
24use mom_domains, only : pass_var,pass_vector
25use mom_error_handler, only : mom_error, fatal, note, warning, is_root_pe
26use mom_file_parser, only : get_param, log_param, log_version, param_file_type
28use mom_grid, only : ocean_grid_type
29use mom_io, only : vardesc, var_desc
31use mom_remapping, only : remappingschemesdoc
32use mom_restart, only : register_restart_field, register_restart_pair, mom_restart_cs
33use mom_restart, only : restart_init, save_restart, query_initialized
35use mom_time_manager, only : time_type
39
40implicit none ; private
41
42#include <MOM_memory.h>
43
44
45! Publicly available functions
49
50! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
51! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
52! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
53! vary with the Boussinesq approximation, the Boussinesq variant is given first.
54
55!> A structure for creating arrays of pointers to 3D arrays with extra gridding information
56type :: p3d
57 integer :: id !< id for FMS external time interpolator
58 integer :: nz_data !< The number of vertical levels in the input field.
59 real, dimension(:,:,:), pointer :: mask_in => null() !< pointer to the data mask (perhaps unused) [nondim]
60 real, dimension(:,:,:), pointer :: p => null() !< pointer to the data, in units that depend
61 !! on the field it refers to [various].
62 real, dimension(:,:,:), pointer :: h => null() !< pointer to the data grid (perhaps unused)
63 !! in [H ~> m or kg m-2]
64end type p3d
65
66!> oda incupd control structure
67type, public :: oda_incupd_cs ; private
68 integer :: nz !< The total number of layers.
69 integer :: nz_data !< The total number of arbritary layers (used by older code).
70 integer :: fldno = 0 !< The number of fields which have already been
71 !! registered by calls to set_up_oda_incupd_field
72
73 type(p3d) :: inc(max_fields_) !< The increments to be applied to the field
74 type(p3d) :: inc_u !< The increments to be applied to the u-velocities, with data in [L T-1 ~> m s-1]
75 type(p3d) :: inc_v !< The increments to be applied to the v-velocities, with data in [L T-1 ~> m s-1]
76 type(p3d) :: ref_h !< Vertical grid on which the increments are provided, with data in [H ~> m or kg m-2]
77
78
79 integer :: nstep_incupd !< number of time step for full update
80 real :: ncount = 0.0 !< increment time step counter [nondim]. This could be an integer
81 !! but a real variable works better with the existing restarts.
82 type(remapping_cs) :: remap_cs !< Remapping parameters and work arrays
83 logical :: incupddataongrid !< True if the incupd data are on the model horizontal grid
84 logical :: uv_inc !< use u and v increments
85
86 ! for diagnostics
87 type(diag_ctrl), pointer :: diag !<structure to regulate output
88 ! diagnostic for inc. fields
89 integer :: id_u_oda_inc = -1 !< diagnostic id for zonal velocity inc.
90 integer :: id_v_oda_inc = -1 !< diagnostic id for meridional velocity inc.
91 integer :: id_h_oda_inc = -1 !< diagnostic id for layer thicknesses inc.
92 integer :: id_t_oda_inc = -1 !< diagnostic id for temperature inc.
93 integer :: id_s_oda_inc = -1 !< diagnostic id for salinity inc.
94
95end type oda_incupd_cs
96
97contains
98
99!> This subroutine defined the control structure of module and register
100!the time counter to full update in restart
101subroutine initialize_oda_incupd_fixed( G, GV, US, CS, restart_CS)
102 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
103 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
104 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
105 type(oda_incupd_cs), pointer :: cs !< A pointer that is set to point to the control
106 !! structure for this module (in/out).
107 type(mom_restart_cs), intent(inout) :: restart_cs !< MOM restart control struct
108
109 if (associated(cs)) then
110 call mom_error(warning, "initialize_oda_incupd_fixed called with an associated "// &
111 "control structure.")
112 return
113 endif
114 allocate(cs)
115
116 ! initialize time counter
117 cs%ncount = 0.0
118 ! register ncount in restart
119 call register_restart_field(cs%ncount, "oda_incupd_ncount", .false., restart_cs,&
120 "Number of inc. update already done", "N/A")
121end subroutine initialize_oda_incupd_fixed
122
123
124!> This subroutine defined the number of time step for full update, stores the layer pressure
125!! increments and initialize remap structure.
126subroutine initialize_oda_incupd( G, GV, US, param_file, CS, data_h, nz_data, restart_CS)
127 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
128 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
129 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
130 integer, intent(in) :: nz_data !< The total number of incr. input layers.
131 type(param_file_type), intent(in) :: param_file !< A structure indicating the open file
132 !! to parse for model parameter values.
133 type(oda_incupd_cs), pointer :: cs !< A pointer that is set to point to the control
134 !! structure for this module (in/out).
135 real, dimension(SZI_(G),SZJ_(G),nz_data), intent(in) :: data_h !< The ODA h
136 !! [H ~> m or kg m-2].
137 type(mom_restart_cs), intent(in) :: restart_cs !< MOM restart control struct
138
139 ! This include declares and sets the variable "version".
140# include "version_variable.h"
141 character(len=40) :: mdl = "MOM_oda" ! This module's name.
142 logical :: use_oda_incupd
143 logical :: bndextrapolation = .true. ! If true, extrapolate boundaries
144 logical :: reset_ncount
145 integer :: i, j, k
146 real :: incupd_timescale ! The amount of timer over which to apply the full update [T ~> s]
147 real :: dt, dt_therm ! Model timesteps [T ~> s]
148 character(len=256) :: mesg
149 character(len=64) :: remapscheme
150 logical :: om4_remap_via_sub_cells ! If true, use the OM4 remapping algorithm
151 real :: h_neglect, h_neglect_edge ! Negligible thicknesses [H ~> m or kg m-2]
152
153 if (.not.associated(cs)) then
154 call mom_error(warning, "initialize_oda_incupd called without an associated "// &
155 "control structure.")
156 return
157 endif
158
159! Set default, read and log parameters
160 call log_version(param_file, mdl, version, "")
161 call get_param(param_file, mdl, "ODA_INCUPD", use_oda_incupd, &
162 "If true, oda incremental updates will be applied "//&
163 "everywhere in the domain.", default=.false.)
164
165 if (.not.use_oda_incupd) return
166
167 call get_param(param_file, mdl, "ODA_INCUPD_NHOURS", incupd_timescale, &
168 "Number of hours for full update (0=direct insertion).", &
169 default=3.0, units="h", scale=3600.0*us%s_to_T)
170 call get_param(param_file, mdl, "ODA_INCUPD_RESET_NCOUNT", reset_ncount, &
171 "If True, reinitialize number of updates already done, ncount.", &
172 default=.true.)
173 call get_param(param_file, mdl, "DT", dt, &
174 "The (baroclinic) dynamics time step. The time-step that "//&
175 "is actually used will be an integer fraction of the "//&
176 "forcing time-step (DT_FORCING in ocean-only mode or the "//&
177 "coupling timestep in coupled mode.)", units="s", scale=us%s_to_T, &
178 fail_if_missing=.true.)
179 call get_param(param_file, mdl, "DT_THERM", dt_therm, &
180 "The thermodynamic and tracer advection time step. "//&
181 "Ideally DT_THERM should be an integer multiple of DT "//&
182 "and less than the forcing or coupling time-step, unless "//&
183 "THERMO_SPANS_COUPLING is true, in which case DT_THERM "//&
184 "can be an integer multiple of the coupling timestep. By "//&
185 "default DT_THERM is set to DT.", &
186 units="s", scale=us%s_to_T, default=us%T_to_s*dt)
187 call get_param(param_file, mdl, "ODA_INCUPD_UV", cs%uv_inc, &
188 "use U,V increments.", &
189 default=.true.)
190 call get_param(param_file, mdl, "REMAPPING_SCHEME", remapscheme, &
191 default="PLM", do_not_log=.true.)
192 call get_param(param_file, mdl, "ODA_REMAPPING_SCHEME", remapscheme, &
193 "This sets the reconstruction scheme used "//&
194 "for vertical remapping for all ODA variables. "//&
195 "It can be one of the following schemes: "//&
196 trim(remappingschemesdoc), default=remapscheme)
197
198 !The default should be REMAP_BOUNDARY_EXTRAP
199 call get_param(param_file, mdl, "BOUNDARY_EXTRAPOLATION", bndextrapolation, &
200 default=.false., do_not_log=.true.)
201 call get_param(param_file, mdl, "ODA_BOUNDARY_EXTRAP", bndextrapolation, &
202 "If true, values at the interfaces of boundary cells are "//&
203 "extrapolated instead of piecewise constant", default=bndextrapolation)
204 call get_param(param_file, mdl, "ODA_INCUPD_DATA_ONGRID", cs%incupdDataOngrid, &
205 "When defined, the incoming oda_incupd data are "//&
206 "assumed to be on the model horizontal grid " , &
207 default=.true.)
208 call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
209 do_not_log=.true., default=.true.)
210 call get_param(param_file, mdl, "ODA_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
211 "If true, use the OM4 remapping-via-subcells algorithm for ODA. "//&
212 "See REMAPPING_USE_OM4_SUBCELLS for more details. "//&
213 "We recommend setting this option to false.", default=om4_remap_via_sub_cells)
214
215 cs%nz = gv%ke
216
217 ! increments on horizontal grid
218 if (.not. cs%incupdDataOngrid) call mom_error(fatal,'initialize_oda_incupd: '// &
219 'The oda_incupd code only applies ODA increments on the same horizontal grid. ')
220
221 ! get number of timestep for full update
222 if (incupd_timescale == 0) then
223 cs%nstep_incupd = 1 !! direct insertion
224 else
225 cs%nstep_incupd = floor( incupd_timescale / dt_therm + 0.001 ) - 1
226 endif
227 write(mesg,'(i12)') cs%nstep_incupd
228 if (is_root_pe()) &
229 call mom_error(note, "initialize_oda_incupd: Number of Timestep of inc. update: "//&
230 trim(mesg))
231
232 ! number of inc. update already done, CS%ncount, either from restart or set to 0.0
233 if (query_initialized(cs%ncount, "oda_incupd_ncount", restart_cs) .and. &
234 .not.reset_ncount) then
235 cs%ncount = cs%ncount
236 else
237 cs%ncount = 0.0
238 endif
239 write(mesg,'(f4.1)') cs%ncount
240 if (is_root_pe()) &
241 call mom_error(note, "initialize_oda_incupd: Inc. update already done: "//&
242 trim(mesg))
243
244 ! get the vertical grid (h_obs) of the increments
245 cs%nz_data = nz_data
246 allocate(cs%Ref_h%p(g%isd:g%ied,g%jsd:g%jed,cs%nz_data), source=0.0)
247 do j=g%jsc,g%jec ; do i=g%isc,g%iec ; do k=1,cs%nz_data
248 cs%Ref_h%p(i,j,k) = data_h(i,j,k)
249 enddo ; enddo ; enddo
250 !### Doing a halo update here on CS%Ref_h%p would avoid needing halo updates each timestep.
251
252 ! Call the constructor for remapping control structure
253 !### Revisit this hard-coded answer_date.
254 if (gv%Boussinesq) then
255 h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
256 else
257 h_neglect = gv%H_subroundoff ; h_neglect_edge = gv%H_subroundoff
258 endif
259
260 call initialize_remapping(cs%remap_cs, remapscheme, boundary_extrapolation=bndextrapolation, &
261 om4_remap_via_sub_cells=om4_remap_via_sub_cells, answer_date=20190101, &
262 h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)
263end subroutine initialize_oda_incupd
264
265
266!> This subroutine stores the increments at h points for the variable
267!! whose address is given by f_ptr.
268subroutine set_up_oda_incupd_field(sp_val, G, GV, CS)
269 type(ocean_grid_type), intent(in) :: g !< Grid structure
270 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
271 type(oda_incupd_cs), pointer :: cs !< oda_incupd control structure (in/out).
272 real, dimension(SZI_(G),SZJ_(G),CS%nz_data), &
273 intent(in) :: sp_val !< increment field, it can have an arbitrary number
274 !! of layers, in various units depending on the
275 !! field it refers to [various].
276
277 integer :: i, j, k
278 character(len=256) :: mesg ! String for error messages
279
280 if (.not.associated(cs)) return
281
282 cs%fldno = cs%fldno + 1
283 if (cs%fldno > max_fields_) then
284 write(mesg,'("Increase MAX_FIELDS_ to at least ",I0," in MOM_memory.h or decrease &
285 &the number of fields increments in the call to &
286 &initialize_oda_incupd." )') cs%fldno
287 call mom_error(fatal,"set_up_oda_incupd_field: "//mesg)
288 endif
289
290 ! store the increment/full field tracer profiles
291 cs%Inc(cs%fldno)%nz_data = cs%nz_data
292 allocate(cs%Inc(cs%fldno)%p(g%isd:g%ied,g%jsd:g%jed,cs%nz_data), source=0.0)
293 do k=1,cs%nz_data ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
294 cs%Inc(cs%fldno)%p(i,j,k) = sp_val(i,j,k)
295 enddo ; enddo ; enddo
296
297end subroutine set_up_oda_incupd_field
298
299
300!> This subroutine stores the increments at u and v points for the variable
301!! whose address is given by u_ptr and v_ptr.
302subroutine set_up_oda_incupd_vel_field(u_val, v_val, G, GV, CS)
303 type(ocean_grid_type), intent(in) :: g !< Grid structure (in).
304 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
305 type(oda_incupd_cs), pointer :: cs !< oda incupd structure (in/out).
306
307 real, dimension(SZIB_(G),SZJ_(G),CS%nz_data), &
308 intent(in) :: u_val !< u increment, it has arbritary number of layers but
309 !! not to exceed the total number of model layers [L T-1 ~> m s-1]
310 real, dimension(SZI_(G),SZJB_(G),CS%nz_data), &
311 intent(in) :: v_val !< v increment, it has arbritary number of layers but
312 !! not to exceed the number of model layers [L T-1 ~> m s-1]
313 integer :: i, j, k
314
315 if (.not.associated(cs)) return
316
317
318 ! store the increment/full field u profile
319 allocate(cs%Inc_u%p(g%isdB:g%iedB,g%jsd:g%jed,cs%nz_data), source=0.0)
320 do j=g%jsc,g%jec ; do i=g%iscB,g%iecB
321 do k=1,cs%nz_data
322 cs%Inc_u%p(i,j,k) = u_val(i,j,k)
323 enddo
324 enddo ; enddo
325
326 ! store the increment/full field v profile
327 allocate(cs%Inc_v%p(g%isd:g%ied,g%jsdB:g%jedB,cs%nz_data), source=0.0)
328 do j=g%jscB,g%jecB ; do i=g%isc,g%iec
329 do k=1,cs%nz_data
330 cs%Inc_v%p(i,j,k) = v_val(i,j,k)
331 enddo
332 enddo ; enddo
333
334end subroutine set_up_oda_incupd_vel_field
335
336! calculation of the increments if using full fields (ODA_INCUPD_INC=.false.) at initialization
337subroutine calc_oda_increments(h, tv, u, v, G, GV, US, CS)
338
339 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure (in).
340 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
341 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
342 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
343 intent(inout) :: h !< Layer thickness [H ~> m or kg m-2] (in)
344 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables
345
346 real, target, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
347 intent(in) :: u !< The zonal velocity that is being
348 !! initialized [L T-1 ~> m s-1]
349 real, target, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
350 intent(in) :: v !< The meridional velocity that is being
351 !! initialized [L T-1 ~> m s-1]
352 type(oda_incupd_cs), pointer :: cs !< A pointer to the control structure for this module
353 !! that is set by a previous call to initialize_oda_incupd (in).
354
355
356 real, dimension(SZK_(GV)) :: tmp_val1 ! data values on the model grid, in rescaled units
357 ! like [S ~> ppt] for salinity.
358 real, allocatable, dimension(:) :: tmp_val2 ! data values remapped to increment grid, in rescaled units
359 ! like [S ~> ppt] for salinity.
360 real, allocatable, dimension(:,:,:) :: h_obs !< Layer-thicknesses of increments [H ~> m or kg m-2]
361 real, allocatable, dimension(:) :: tmp_h ! temporary array for corrected h_obs [H ~> m or kg m-2]
362 real, allocatable, dimension(:) :: hu_obs ! A column of observation-grid thicknesses at u points [H ~> m or kg m-2]
363 real, allocatable, dimension(:) :: hv_obs ! A column of observation-grid thicknesses at v points [H ~> m or kg m-2]
364 real, dimension(SZK_(GV)) :: hu, hv ! A column of thicknesses at u or v points [H ~> m or kg m-2]
365
366
367 integer :: i, j, k, is, ie, js, je, nz, nz_data
368 integer :: isb, ieb, jsb, jeb
369 real :: sum_h1, sum_h2 ! vertical sums of h's [H ~> m or kg m-2]
370
371 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
372 isb = g%iscB ; ieb = g%iecB ; jsb = g%jscB ; jeb = g%jecB
373 if (.not.associated(cs)) return
374
375
376 ! increments calculated on if CS%ncount = 0.0
377 if (cs%ncount /= 0.0) call mom_error(fatal,'calc_oda_increments: '// &
378 'CS%ncount should be 0.0 to get accurate increments.')
379
380 ! get h_obs
381 nz_data = cs%Inc(1)%nz_data
382 allocate(h_obs(g%isd:g%ied,g%jsd:g%jed,nz_data), source=0.0)
383 do k=1,nz_data ; do j=js,je ; do i=is,ie
384 h_obs(i,j,k) = cs%Ref_h%p(i,j,k)
385 enddo ; enddo ; enddo
386 call pass_var(h_obs,g%Domain)
387
388
389 ! allocate 1-d arrays
390 allocate(tmp_h(nz_data), source=0.0)
391 allocate(tmp_val2(nz_data), source=0.0)
392 allocate(hu_obs(nz_data), source=0.0)
393 allocate(hv_obs(nz_data), source=0.0)
394
395 ! remap t,s (on h_init) to h_obs to get increment
396 tmp_val1(:) = 0.0
397 do j=js,je ; do i=is,ie
398 if (g%mask2dT(i,j) == 1) then
399 ! account for the different SSH
400 sum_h1 = 0.0
401 sum_h2 = 0.0
402 do k=1,nz
403 sum_h1 = sum_h1+h(i,j,k)
404 enddo
405
406 do k=1,nz_data
407 sum_h2 = sum_h2+h_obs(i,j,k)
408 enddo
409 do k=1,nz_data
410 tmp_h(k)=(sum_h1/sum_h2)*h_obs(i,j,k)
411 enddo
412 ! get temperature
413 do k=1,nz
414 tmp_val1(k) = tv%T(i,j,k)
415 enddo
416 ! remap tracer on h_obs
417 call remapping_core_h(cs%remap_cs, nz, h(i,j,1:nz), tmp_val1, &
418 nz_data, tmp_h(1:nz_data), tmp_val2)
419 ! get increment from full field on h_obs
420 do k=1,nz_data
421 cs%Inc(1)%p(i,j,k) = cs%Inc(1)%p(i,j,k) - tmp_val2(k)
422 enddo
423
424 ! get salinity
425 do k=1,nz
426 tmp_val1(k) = tv%S(i,j,k)
427 enddo
428 ! remap tracer on h_obs
429 call remapping_core_h(cs%remap_cs, nz, h(i,j,1:nz), tmp_val1, &
430 nz_data, tmp_h(1:nz_data), tmp_val2)
431 ! get increment from full field on h_obs
432 do k=1,nz_data
433 cs%Inc(2)%p(i,j,k) = cs%Inc(2)%p(i,j,k) - tmp_val2(k)
434 enddo
435 endif
436 enddo ; enddo
437
438 ! remap u to h_obs to get increment
439 if (cs%uv_inc) then
440 call pass_var(h, g%Domain)
441
442 hu(:) = 0.0
443 do j=js,je ; do i=isb,ieb
444 if (g%mask2dCu(i,j) == 1) then
445 ! get u-velocity
446 do k=1,nz
447 tmp_val1(k) = u(i,j,k)
448 ! get the h and h_obs at u points
449 hu(k) = 0.5*( h(i,j,k)+ h(i+1,j,k))
450 enddo
451 do k=1,nz_data
452 hu_obs(k) = 0.5*(h_obs(i,j,k)+h_obs(i+1,j,k))
453 enddo
454 ! account for the different SSH
455 sum_h1 = 0.0
456 do k=1,nz
457 sum_h1 = sum_h1+hu(k)
458 enddo
459 sum_h2 = 0.0
460 do k=1,nz_data
461 sum_h2 = sum_h2+hu_obs(k)
462 enddo
463 do k=1,nz_data
464 hu_obs(k)=(sum_h1/sum_h2)*hu_obs(k)
465 enddo
466 ! remap model u on hu_obs
467 call remapping_core_h(cs%remap_cs, nz, hu(1:nz), tmp_val1, &
468 nz_data, hu_obs(1:nz_data), tmp_val2)
469 ! get increment from full field on h_obs
470 do k=1,nz_data
471 cs%Inc_u%p(i,j,k) = cs%Inc_u%p(i,j,k) - tmp_val2(k)
472 enddo
473 endif
474 enddo ; enddo
475
476 ! remap v to h_obs to get increment
477 hv(:) = 0.0
478 do j=jsb,jeb ; do i=is,ie
479 if (g%mask2dCv(i,j) == 1) then
480 ! get v-velocity
481 do k=1,nz
482 tmp_val1(k) = v(i,j,k)
483 ! get the h and h_obs at v points
484 hv(k) = 0.5*(h(i,j,k)+h(i,j+1,k))
485 enddo
486 do k=1,nz_data
487 hv_obs(k) = 0.5*(h_obs(i,j,k)+h_obs(i,j+1,k))
488 enddo
489 ! account for the different SSH
490 sum_h1 = 0.0
491 do k=1,nz
492 sum_h1 = sum_h1+hv(k)
493 enddo
494 sum_h2 = 0.0
495 do k=1,nz_data
496 sum_h2 = sum_h2+hv_obs(k)
497 enddo
498 do k=1,nz_data
499 hv_obs(k)=(sum_h1/sum_h2)*hv_obs(k)
500 enddo
501 ! remap model v on hv_obs
502 call remapping_core_h(cs%remap_cs, nz, hv(1:nz), tmp_val1, &
503 nz_data, hv_obs(1:nz_data), tmp_val2)
504 ! get increment from full field on h_obs
505 do k=1,nz_data
506 cs%Inc_v%p(i,j,k) = cs%Inc_v%p(i,j,k) - tmp_val2(k)
507 enddo
508 endif
509 enddo ; enddo
510 endif ! uv_inc
511
512 call pass_var(cs%Inc(1)%p, g%Domain)
513 call pass_var(cs%Inc(2)%p, g%Domain)
514 call pass_vector(cs%Inc_u%p,cs%Inc_v%p,g%Domain)
515
516 ! deallocate arrays
517 deallocate(tmp_h,tmp_val2,hu_obs,hv_obs)
518 deallocate(h_obs)
519
520end subroutine calc_oda_increments
521
522!> This subroutine applies oda increments to layers thicknesses, temp, salt, U
523!! and V everywhere .
524subroutine apply_oda_incupd(h, tv, u, v, dt, G, GV, US, CS)
525 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure (in).
526 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
527 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
528 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
529 intent(inout) :: h !< Layer thickness [H ~> m or kg m-2] (in)
530 type(thermo_var_ptrs), intent(inout) :: tv !< A structure pointing to various thermodynamic variables
531
532 real, target, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
533 intent(inout) :: u !< The zonal velocity that is being
534 !! initialized [L T-1 ~> m s-1]
535 real, target, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
536 intent(inout) :: v !< The meridional velocity that is being
537 !! initialized [L T-1 ~> m s-1]
538
539 real, intent(in) :: dt !< The amount of time covered by this call [T ~> s].
540 type(oda_incupd_cs), pointer :: cs !< A pointer to the control structure for this module
541 !! that is set by a previous call to initialize_oda_incupd (in).
542
543 ! Local variables
544 real, allocatable, dimension(:) :: tmp_val2 ! data values remapped to increment grid, in rescaled units
545 ! like [S ~> ppt] for salinity.
546 real, dimension(SZK_(GV)) :: tmp_val1 ! data values on the model grid, in rescaled units
547 ! like [S ~> ppt] for salinity.
548 real, dimension(SZK_(GV)) :: hu, hv ! A column of thicknesses at u or v points [H ~> m or kg m-2]
549
550 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: tmp_t !< A temporary array for t increments [C ~> degC]
551 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: tmp_s !< A temporary array for s increments [S ~> ppt]
552 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: tmp_u !< A temporary array for u increments [L T-1 ~> m s-1]
553 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: tmp_v !< A temporary array for v increments [L T-1 ~> m s-1]
554
555 real, allocatable, dimension(:,:,:) :: h_obs !< h of increments [H ~> m or kg m-2]
556 real, allocatable, dimension(:) :: tmp_h !< temporary array for corrected h_obs [H ~> m or kg m-2]
557 real, allocatable, dimension(:) :: hu_obs ! A column of observation-grid thicknesses at u points [H ~> m or kg m-2]
558 real, allocatable, dimension(:) :: hv_obs ! A column of observation-grid thicknesses at v points [H ~> m or kg m-2]
559
560 integer :: i, j, k, is, ie, js, je, nz, nz_data
561 integer :: isb, ieb, jsb, jeb
562! integer :: ncount ! time step counter
563 real :: inc_wt ! weight of the update for this time-step [nondim]
564 real :: sum_h1, sum_h2 ! vertical sums of h's [H ~> m or kg m-2]
565 character(len=256) :: mesg
566
567 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
568 isb = g%iscB ; ieb = g%iecB ; jsb = g%jscB ; jeb = g%jecB
569 if (.not.associated(cs)) return
570
571 ! no assimilation after CS%step_incupd
572 if (cs%ncount >= cs%nstep_incupd) then
573 if (is_root_pe()) call mom_error(note,"ended updating fields with increments. ")
574 return
575 endif !ncount>CS%nstep_incupd
576
577 ! update counter
578 cs%ncount = cs%ncount+1.0
579 inc_wt = 1.0/cs%nstep_incupd
580
581 ! print out increments
582 write(mesg,'(f10.0)') cs%ncount
583 if (is_root_pe()) call mom_error(note,"updating fields with increments ncount:"//trim(mesg))
584 write(mesg,'(f10.8)') inc_wt
585 if (is_root_pe()) call mom_error(note,"updating fields with weight inc_wt:"//trim(mesg))
586
587 ! get h_obs
588 nz_data = cs%Inc(1)%nz_data
589 allocate(h_obs(g%isd:g%ied,g%jsd:g%jed,nz_data), source=0.0)
590 do k=1,nz_data ; do j=js,je ; do i=is,ie
591 h_obs(i,j,k) = cs%Ref_h%p(i,j,k)
592 enddo ; enddo ; enddo
593 call pass_var(h_obs,g%Domain)
594
595 ! allocate 1-d array
596 allocate(tmp_h(nz_data), source=0.0)
597 allocate(tmp_val2(nz_data))
598 allocate(hu_obs(nz_data), source=0.0)
599 allocate(hv_obs(nz_data), source=0.0)
600
601 ! add increments to tracers
602 tmp_val1(:) = 0.0
603 tmp_t(:,:,:) = 0.0 ; tmp_s(:,:,:) = 0.0 ! diagnostics
604 do j=js,je ; do i=is,ie
605 if (g%mask2dT(i,j) == 1) then
606 ! account for the different SSH
607 sum_h1 = 0.0
608 do k=1,nz
609 sum_h1 = sum_h1+h(i,j,k)
610 enddo
611 sum_h2 = 0.0
612 do k=1,nz_data
613 sum_h2 = sum_h2+h_obs(i,j,k)
614 enddo
615 do k=1,nz_data
616 tmp_h(k) = ( sum_h1 / sum_h2 ) * h_obs(i,j,k)
617 enddo
618 ! get temperature increment
619 do k=1,nz_data
620 tmp_val2(k) = cs%Inc(1)%p(i,j,k)
621 enddo
622 ! remap increment profile on model h
623 call remapping_core_h(cs%remap_cs, nz_data, tmp_h(1:nz_data), tmp_val2, &
624 nz, h(i,j,1:nz), tmp_val1)
625 do k=1,nz
626 ! add increment to tracer on model h
627 tv%T(i,j,k) = tv%T(i,j,k) + inc_wt * tmp_val1(k)
628 tmp_t(i,j,k) = tmp_val1(k) ! store T increment for diagnostics
629 enddo
630
631 ! get salinity increment
632 do k=1,nz_data
633 tmp_val2(k) = cs%Inc(2)%p(i,j,k)
634 enddo
635 ! remap increment profile on model h
636 call remapping_core_h(cs%remap_cs, nz_data, tmp_h(1:nz_data), tmp_val2, &
637 nz, h(i,j,1:nz), tmp_val1)
638 ! add increment to tracer on model h
639 do k=1,nz
640 tv%S(i,j,k) = tv%S(i,j,k) + inc_wt * tmp_val1(k)
641 tmp_s(i,j,k) = tmp_val1(k) ! store S increment for diagnostics
642 ! bound salinity values ! check if it is correct to do that or if it hides
643 ! other problems ...
644 tv%S(i,j,k) = max(0.0 , tv%S(i,j,k))
645 enddo
646 endif
647 enddo ; enddo
648
649
650 ! add u and v increments
651 if (cs%uv_inc) then
652
653 call pass_var(h,g%Domain) ! to ensure reproducibility
654
655 ! add increments to u
656 hu(:) = 0.0
657 tmp_u(:,:,:) = 0.0 ! diagnostics
658 do j=js,je ; do i=isb,ieb
659 if (g%mask2dCu(i,j) == 1) then
660 do k=1,nz_data
661 ! get u increment
662 tmp_val2(k) = cs%Inc_u%p(i,j,k)
663 ! get the h and h_obs at u points
664 hu_obs(k) = 0.5 * ( h_obs(i,j,k) + h_obs(i+1,j,k) )
665 enddo
666 do k=1,nz
667 hu(k) = 0.5 * ( h(i,j,k) + h(i+1,j,k) )
668 enddo
669 ! account for different SSH
670 sum_h1 = 0.0
671 do k=1,nz
672 sum_h1 = sum_h1 + hu(k)
673 enddo
674 sum_h2 = 0.0
675 do k=1,nz_data
676 sum_h2 = sum_h2 + hu_obs(k)
677 enddo
678 do k=1,nz_data
679 hu_obs(k) = ( sum_h1 / sum_h2 ) * hu_obs(k)
680 enddo
681 ! remap increment profile on hu
682 call remapping_core_h(cs%remap_cs, nz_data, hu_obs(1:nz_data), tmp_val2, &
683 nz, hu(1:nz), tmp_val1)
684 ! add increment to u-velocity on hu
685 do k=1,nz
686 u(i,j,k) = u(i,j,k) + inc_wt * tmp_val1(k)
687 ! store increment for diagnostics
688 tmp_u(i,j,k) = tmp_val1(k)
689 enddo
690 endif
691 enddo ; enddo
692
693 ! add increments to v
694 hv(:) = 0.0
695 tmp_v(:,:,:) = 0.0 ! diagnostics
696 do j=jsb,jeb ; do i=is,ie
697 if (g%mask2dCv(i,j) == 1) then
698 ! get v increment
699 do k=1,nz_data
700 tmp_val2(k) = cs%Inc_v%p(i,j,k)
701 ! get the h and h_obs at v points
702 hv_obs(k) = 0.5 * ( h_obs(i,j,k) + h_obs(i,j+1,k) )
703 enddo
704 do k=1,nz
705 hv(k) = 0.5 * (h(i,j,k) + h(i,j+1,k) )
706 enddo
707 ! account for different SSH
708 sum_h1 = 0.0
709 do k=1,nz
710 sum_h1 = sum_h1 + hv(k)
711 enddo
712 sum_h2 = 0.0
713 do k=1,nz_data
714 sum_h2 = sum_h2 + hv_obs(k)
715 enddo
716 do k=1,nz_data
717 hv_obs(k)=( sum_h1 / sum_h2 ) * hv_obs(k)
718 enddo
719 ! remap increment profile on hv
720 call remapping_core_h(cs%remap_cs, nz_data, hv_obs(1:nz_data), tmp_val2, &
721 nz, hv(1:nz), tmp_val1)
722 ! add increment to v-velocity on hv
723 do k=1,nz
724 v(i,j,k) = v(i,j,k) + inc_wt * tmp_val1(k)
725 ! store increment for diagnostics
726 tmp_v(i,j,k) = tmp_val1(k)
727 enddo
728 endif
729 enddo ; enddo
730
731 endif ! uv_inc
732
733 call pass_var(tv%T, g%Domain)
734 call pass_var(tv%S, g%Domain)
735 call pass_vector(u,v,g%Domain)
736
737 ! Diagnostics of increments, mostly used for debugging.
738 if (cs%uv_inc) then
739 if (cs%id_u_oda_inc > 0) call post_data(cs%id_u_oda_inc, tmp_u, cs%diag)
740 if (cs%id_v_oda_inc > 0) call post_data(cs%id_v_oda_inc, tmp_v, cs%diag)
741 endif
742 !### The argument here seems wrong.
743 if (cs%id_h_oda_inc > 0) call post_data(cs%id_h_oda_inc, h , cs%diag)
744 if (cs%id_T_oda_inc > 0) call post_data(cs%id_T_oda_inc, tmp_t, cs%diag)
745 if (cs%id_S_oda_inc > 0) call post_data(cs%id_S_oda_inc, tmp_s, cs%diag)
746
747 ! deallocate arrays
748 deallocate(tmp_h,tmp_val2,hu_obs,hv_obs)
749 deallocate(h_obs)
750
751end subroutine apply_oda_incupd
752
753!> Output increment if using full fields for the oda_incupd module.
754subroutine output_oda_incupd_inc(Time, G, GV, param_file, CS, US)
755 type(time_type), target, intent(in) :: time !< The current model time
756 type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure
757 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
758 type(param_file_type), intent(in) :: param_file !< A structure indicating the open file
759 !! to parse for
760 !model parameter
761 !values.
762 type(oda_incupd_cs), pointer :: cs !< ODA incupd control structure
763 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling
764
765 type(mom_restart_cs), pointer :: restart_csp_tmp => null()
766
767 type(directories) :: dirs
768 type(vardesc) :: u_desc, v_desc
769
770 character(len=40) :: mdl = "MOM_oda" ! This module's name.
771 character(len=200) :: inc_file ! name of the increment file
772
773 if (.not.associated(cs)) return
774 ! get the output_directory
775 call get_mom_input(dirs=dirs)
776 if (is_root_pe()) call mom_error(note,"output increments in output_directory")
777
778 ! get a restart structure
779 call restart_init(param_file, restart_csp_tmp)
780
781 ! register the variables to write
782 call register_restart_field(cs%Inc(1)%p, "T_inc", .true., restart_csp_tmp, &
783 "Pot. T. increment", "degC", conversion=us%C_to_degC)
784 call register_restart_field(cs%Inc(2)%p, "S_inc", .true., restart_csp_tmp, &
785 "Salinity increment", "psu", conversion=us%S_to_ppt)
786 call register_restart_field(cs%Ref_h%p, "h_obs", .true., restart_csp_tmp, &
787 "Observational h", units=get_thickness_units(gv), conversion=gv%H_to_MKS)
788 if (cs%uv_inc) then
789 u_desc = var_desc("u_inc", "m s-1", "U-vel increment", hor_grid='Cu')
790 v_desc = var_desc("v_inc", "m s-1", "V-vel increment", hor_grid='Cv')
791 call register_restart_pair(cs%Inc_u%p, cs%Inc_v%p, u_desc, v_desc, &
792 .false., restart_csp_tmp, conversion=us%L_T_to_m_s)
793 endif
794
795 ! get the name of the output file
796 call get_param(param_file, mdl, "ODA_INCUPD_OUTPUT_FILE", inc_file,&
797 "The name-root of the output file for the increment if using full fields.", &
798 default="MOM.inc")
799
800 ! write the increments file
801 call save_restart(dirs%output_directory, time, g, restart_csp_tmp, &
802 filename=inc_file, gv=gv) !, write_ic=.true.)
803
804end subroutine output_oda_incupd_inc
805
806
807
808!> Initialize diagnostics for the oda_incupd module.
809subroutine init_oda_incupd_diags(Time, G, GV, diag, CS, US)
810 type(time_type), target, intent(in) :: time !< The current model time
811 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
812 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
813 type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic
814 !! output.
815 type(oda_incupd_cs), pointer :: cs !< ALE sponge control structure
816 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling
817
818 if (.not.associated(cs)) return
819
820 cs%diag => diag
821 ! These diagnostics of the state variables increments are useful for debugging the ODA code.
822 cs%id_u_oda_inc = register_diag_field('ocean_model', 'u_oda_inc', diag%axesCuL, time, &
823 'Zonal velocity ODA inc.', 'm s-1', conversion=us%L_T_to_m_s)
824 cs%id_v_oda_inc = register_diag_field('ocean_model', 'v_oda_inc', diag%axesCvL, time, &
825 'Meridional velocity ODA inc.', 'm s-1', conversion=us%L_T_to_m_s)
826 cs%id_h_oda_inc = register_diag_field('ocean_model', 'h_oda_inc', diag%axesTL, time, &
827 'Layer Thickness ODA inc.', get_thickness_units(gv), conversion=gv%H_to_mks)
828 cs%id_T_oda_inc = register_diag_field('ocean_model', 'T_oda_inc', diag%axesTL, time, &
829 'Temperature ODA inc.', 'degC', conversion=us%C_to_degC)
830 cs%id_S_oda_inc = register_diag_field('ocean_model', 'S_oda_inc', diag%axesTL, time, &
831 'Salinity ODA inc.', 'PSU', conversion=us%S_to_ppt)
832
833end subroutine init_oda_incupd_diags
834
835!> This subroutine deallocates any memory associated with the oda_incupd module.
836subroutine oda_incupd_end(CS)
837 type(oda_incupd_cs), pointer :: cs !< A pointer to the control structure that is
838 !! set by a previous call to initialize_oda_incupd.
839
840 integer :: m
841
842 if (.not.associated(cs)) return
843
844 do m=1,cs%fldno
845 if (associated(cs%Inc(m)%p)) deallocate(cs%Inc(m)%p)
846 enddo
847
848 deallocate(cs)
849
850end subroutine oda_incupd_end
851
852end module mom_oda_incupd