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%nstep_incupd
572 if (cs%ncount >= cs%nstep_incupd) then
573 return
574 endif !ncount>CS%nstep_incupd
575
576 ! update counter
577 cs%ncount = cs%ncount+1.0
578 inc_wt = 1.0/cs%nstep_incupd
579
580 ! print out increments
581 write(mesg,'(f10.0)') cs%ncount
582 if (is_root_pe()) call mom_error(note,"updating fields with increments ncount:"//trim(mesg))
583 write(mesg,'(f10.8)') inc_wt
584 if (is_root_pe()) call mom_error(note,"updating fields with weight inc_wt:"//trim(mesg))
585
586 ! get h_obs
587 nz_data = cs%Inc(1)%nz_data
588 allocate(h_obs(g%isd:g%ied,g%jsd:g%jed,nz_data), source=0.0)
589 do k=1,nz_data ; do j=js,je ; do i=is,ie
590 h_obs(i,j,k) = cs%Ref_h%p(i,j,k)
591 enddo ; enddo ; enddo
592 call pass_var(h_obs,g%Domain)
593
594 ! allocate 1-d array
595 allocate(tmp_h(nz_data), source=0.0)
596 allocate(tmp_val2(nz_data))
597 allocate(hu_obs(nz_data), source=0.0)
598 allocate(hv_obs(nz_data), source=0.0)
599
600 ! add increments to tracers
601 tmp_val1(:) = 0.0
602 tmp_t(:,:,:) = 0.0 ; tmp_s(:,:,:) = 0.0 ! diagnostics
603 do j=js,je ; do i=is,ie
604 if (g%mask2dT(i,j) == 1) then
605 ! account for the different SSH
606 sum_h1 = 0.0
607 do k=1,nz
608 sum_h1 = sum_h1+h(i,j,k)
609 enddo
610 sum_h2 = 0.0
611 do k=1,nz_data
612 sum_h2 = sum_h2+h_obs(i,j,k)
613 enddo
614 do k=1,nz_data
615 tmp_h(k) = ( sum_h1 / sum_h2 ) * h_obs(i,j,k)
616 enddo
617 ! get temperature increment
618 do k=1,nz_data
619 tmp_val2(k) = cs%Inc(1)%p(i,j,k)
620 enddo
621 ! remap increment profile on model h
622 call remapping_core_h(cs%remap_cs, nz_data, tmp_h(1:nz_data), tmp_val2, &
623 nz, h(i,j,1:nz), tmp_val1)
624 do k=1,nz
625 ! add increment to tracer on model h
626 tv%T(i,j,k) = tv%T(i,j,k) + inc_wt * tmp_val1(k)
627 tmp_t(i,j,k) = tmp_val1(k) ! store T increment for diagnostics
628 enddo
629
630 ! get salinity increment
631 do k=1,nz_data
632 tmp_val2(k) = cs%Inc(2)%p(i,j,k)
633 enddo
634 ! remap increment profile on model h
635 call remapping_core_h(cs%remap_cs, nz_data, tmp_h(1:nz_data), tmp_val2, &
636 nz, h(i,j,1:nz), tmp_val1)
637 ! add increment to tracer on model h
638 do k=1,nz
639 tv%S(i,j,k) = tv%S(i,j,k) + inc_wt * tmp_val1(k)
640 tmp_s(i,j,k) = tmp_val1(k) ! store S increment for diagnostics
641 ! bound salinity values ! check if it is correct to do that or if it hides
642 ! other problems ...
643 tv%S(i,j,k) = max(0.0 , tv%S(i,j,k))
644 enddo
645 endif
646 enddo ; enddo
647
648
649 ! add u and v increments
650 if (cs%uv_inc) then
651
652 call pass_var(h,g%Domain) ! to ensure reproducibility
653
654 ! add increments to u
655 hu(:) = 0.0
656 tmp_u(:,:,:) = 0.0 ! diagnostics
657 do j=js,je ; do i=isb,ieb
658 if (g%mask2dCu(i,j) == 1) then
659 do k=1,nz_data
660 ! get u increment
661 tmp_val2(k) = cs%Inc_u%p(i,j,k)
662 ! get the h and h_obs at u points
663 hu_obs(k) = 0.5 * ( h_obs(i,j,k) + h_obs(i+1,j,k) )
664 enddo
665 do k=1,nz
666 hu(k) = 0.5 * ( h(i,j,k) + h(i+1,j,k) )
667 enddo
668 ! account for different SSH
669 sum_h1 = 0.0
670 do k=1,nz
671 sum_h1 = sum_h1 + hu(k)
672 enddo
673 sum_h2 = 0.0
674 do k=1,nz_data
675 sum_h2 = sum_h2 + hu_obs(k)
676 enddo
677 do k=1,nz_data
678 hu_obs(k) = ( sum_h1 / sum_h2 ) * hu_obs(k)
679 enddo
680 ! remap increment profile on hu
681 call remapping_core_h(cs%remap_cs, nz_data, hu_obs(1:nz_data), tmp_val2, &
682 nz, hu(1:nz), tmp_val1)
683 ! add increment to u-velocity on hu
684 do k=1,nz
685 u(i,j,k) = u(i,j,k) + inc_wt * tmp_val1(k)
686 ! store increment for diagnostics
687 tmp_u(i,j,k) = tmp_val1(k)
688 enddo
689 endif
690 enddo ; enddo
691
692 ! add increments to v
693 hv(:) = 0.0
694 tmp_v(:,:,:) = 0.0 ! diagnostics
695 do j=jsb,jeb ; do i=is,ie
696 if (g%mask2dCv(i,j) == 1) then
697 ! get v increment
698 do k=1,nz_data
699 tmp_val2(k) = cs%Inc_v%p(i,j,k)
700 ! get the h and h_obs at v points
701 hv_obs(k) = 0.5 * ( h_obs(i,j,k) + h_obs(i,j+1,k) )
702 enddo
703 do k=1,nz
704 hv(k) = 0.5 * (h(i,j,k) + h(i,j+1,k) )
705 enddo
706 ! account for different SSH
707 sum_h1 = 0.0
708 do k=1,nz
709 sum_h1 = sum_h1 + hv(k)
710 enddo
711 sum_h2 = 0.0
712 do k=1,nz_data
713 sum_h2 = sum_h2 + hv_obs(k)
714 enddo
715 do k=1,nz_data
716 hv_obs(k)=( sum_h1 / sum_h2 ) * hv_obs(k)
717 enddo
718 ! remap increment profile on hv
719 call remapping_core_h(cs%remap_cs, nz_data, hv_obs(1:nz_data), tmp_val2, &
720 nz, hv(1:nz), tmp_val1)
721 ! add increment to v-velocity on hv
722 do k=1,nz
723 v(i,j,k) = v(i,j,k) + inc_wt * tmp_val1(k)
724 ! store increment for diagnostics
725 tmp_v(i,j,k) = tmp_val1(k)
726 enddo
727 endif
728 enddo ; enddo
729
730 endif ! uv_inc
731
732 call pass_var(tv%T, g%Domain)
733 call pass_var(tv%S, g%Domain)
734 call pass_vector(u,v,g%Domain)
735
736 ! Diagnostics of increments, mostly used for debugging.
737 if (cs%uv_inc) then
738 if (cs%id_u_oda_inc > 0) call post_data(cs%id_u_oda_inc, tmp_u, cs%diag)
739 if (cs%id_v_oda_inc > 0) call post_data(cs%id_v_oda_inc, tmp_v, cs%diag)
740 endif
741 !### The argument here seems wrong.
742 if (cs%id_h_oda_inc > 0) call post_data(cs%id_h_oda_inc, h , cs%diag)
743 if (cs%id_T_oda_inc > 0) call post_data(cs%id_T_oda_inc, tmp_t, cs%diag)
744 if (cs%id_S_oda_inc > 0) call post_data(cs%id_S_oda_inc, tmp_s, cs%diag)
745
746 ! deallocate arrays
747 deallocate(tmp_h,tmp_val2,hu_obs,hv_obs)
748 deallocate(h_obs)
749
750end subroutine apply_oda_incupd
751
752!> Output increment if using full fields for the oda_incupd module.
753subroutine output_oda_incupd_inc(Time, G, GV, param_file, CS, US)
754 type(time_type), target, intent(in) :: time !< The current model time
755 type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure
756 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
757 type(param_file_type), intent(in) :: param_file !< A structure indicating the open file
758 !! to parse for
759 !model parameter
760 !values.
761 type(oda_incupd_cs), pointer :: cs !< ODA incupd control structure
762 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling
763
764 type(mom_restart_cs), pointer :: restart_csp_tmp => null()
765
766 type(directories) :: dirs
767 type(vardesc) :: u_desc, v_desc
768
769 character(len=40) :: mdl = "MOM_oda" ! This module's name.
770 character(len=200) :: inc_file ! name of the increment file
771
772 if (.not.associated(cs)) return
773 ! get the output_directory
774 call get_mom_input(dirs=dirs)
775 if (is_root_pe()) call mom_error(note,"output increments in output_directory")
776
777 ! get a restart structure
778 call restart_init(param_file, restart_csp_tmp)
779
780 ! register the variables to write
781 call register_restart_field(cs%Inc(1)%p, "T_inc", .true., restart_csp_tmp, &
782 "Pot. T. increment", "degC", conversion=us%C_to_degC)
783 call register_restart_field(cs%Inc(2)%p, "S_inc", .true., restart_csp_tmp, &
784 "Salinity increment", "psu", conversion=us%S_to_ppt)
785 call register_restart_field(cs%Ref_h%p, "h_obs", .true., restart_csp_tmp, &
786 "Observational h", units=get_thickness_units(gv), conversion=gv%H_to_MKS)
787 if (cs%uv_inc) then
788 u_desc = var_desc("u_inc", "m s-1", "U-vel increment", hor_grid='Cu')
789 v_desc = var_desc("v_inc", "m s-1", "V-vel increment", hor_grid='Cv')
790 call register_restart_pair(cs%Inc_u%p, cs%Inc_v%p, u_desc, v_desc, &
791 .false., restart_csp_tmp, conversion=us%L_T_to_m_s)
792 endif
793
794 ! get the name of the output file
795 call get_param(param_file, mdl, "ODA_INCUPD_OUTPUT_FILE", inc_file,&
796 "The name-root of the output file for the increment if using full fields.", &
797 default="MOM.inc")
798
799 ! write the increments file
800 call save_restart(dirs%output_directory, time, g, restart_csp_tmp, &
801 filename=inc_file, gv=gv) !, write_ic=.true.)
802
803end subroutine output_oda_incupd_inc
804
805
806
807!> Initialize diagnostics for the oda_incupd module.
808subroutine init_oda_incupd_diags(Time, G, GV, diag, CS, US)
809 type(time_type), target, intent(in) :: time !< The current model time
810 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
811 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
812 type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic
813 !! output.
814 type(oda_incupd_cs), pointer :: cs !< ALE sponge control structure
815 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling
816
817 if (.not.associated(cs)) return
818
819 cs%diag => diag
820 ! These diagnostics of the state variables increments are useful for debugging the ODA code.
821 cs%id_u_oda_inc = register_diag_field('ocean_model', 'u_oda_inc', diag%axesCuL, time, &
822 'Zonal velocity ODA inc.', 'm s-1', conversion=us%L_T_to_m_s)
823 cs%id_v_oda_inc = register_diag_field('ocean_model', 'v_oda_inc', diag%axesCvL, time, &
824 'Meridional velocity ODA inc.', 'm s-1', conversion=us%L_T_to_m_s)
825 cs%id_h_oda_inc = register_diag_field('ocean_model', 'h_oda_inc', diag%axesTL, time, &
826 'Layer Thickness ODA inc.', get_thickness_units(gv), conversion=gv%H_to_mks)
827 cs%id_T_oda_inc = register_diag_field('ocean_model', 'T_oda_inc', diag%axesTL, time, &
828 'Temperature ODA inc.', 'degC', conversion=us%C_to_degC)
829 cs%id_S_oda_inc = register_diag_field('ocean_model', 'S_oda_inc', diag%axesTL, time, &
830 'Salinity ODA inc.', 'PSU', conversion=us%S_to_ppt)
831
832end subroutine init_oda_incupd_diags
833
834!> This subroutine deallocates any memory associated with the oda_incupd module.
835subroutine oda_incupd_end(CS)
836 type(oda_incupd_cs), pointer :: cs !< A pointer to the control structure that is
837 !! set by a previous call to initialize_oda_incupd.
838
839 integer :: m
840
841 if (.not.associated(cs)) return
842
843 do m=1,cs%fldno
844 if (associated(cs%Inc(m)%p)) deallocate(cs%Inc(m)%p)
845 enddo
846
847 deallocate(cs)
848
849end subroutine oda_incupd_end
850
851end module mom_oda_incupd