MOM_oda_driver.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!> Interfaces for MOM6 ensembles and data assimilation.
7
8! This file is part of MOM6. see LICENSE.md for the license.
9
10! MOM infrastructure
11use mom_coms, only : pe_here, num_pes
12use mom_coms, only : set_pelist, set_rootpe, get_pelist, broadcast
13use mom_domains, only : domain2d, global_field, get_domain_extent
14use mom_domains, only : pass_var, redistribute_array, broadcast_domain
15use mom_diag_mediator, only : register_diag_field, diag_axis_init, post_data
18use mom_ensemble_manager, only : get_ensemble_id, get_ensemble_size
19use mom_ensemble_manager, only : get_ensemble_pelist, get_ensemble_filter_pelist
20use mom_error_handler, only : stdout, stdlog, mom_error
21use mom_io, only : single_file
22use mom_interp_infra, only : init_extern_field
23use mom_interp_infra, only : time_interp_extern
24use mom_interpolate, only : external_field
26use mom_remapping, only : remappingschemesdoc
27use mom_time_manager, only : time_type, real_to_time, get_date
28use mom_time_manager, only : operator(+), operator(>=), operator(/=)
29use mom_time_manager, only : operator(==), operator(<)
30use mom_cpu_clock, only : cpu_clock_begin, cpu_clock_end, cpu_clock_id
31use mom_horizontal_regridding, only : horiz_interp_and_extrap_tracer
32! ODA Modules
35!This preprocessing directive enables the SPEAR online ensemble data assimilation
36!configuration. Existing community based APIs for data assimilation are currently
37!called offline for forecast applications using information read from a MOM6 state file.
38!The SPEAR configuration (https://doi.org/10.1029/2020MS002149) calculated increments
39!efficiently online. A community-based set of APIs should be implemented in place
40!of the CPP directive when this is available.
41#ifdef ENABLE_ECDA
42use eakf_oda_mod, only : ensemble_filter
43#endif
44use kdtree, only : kd_root !# A kd-tree object using JEDI APIs
45! MOM Modules
46use mom_io, only : slasher, mom_read_data
48use mom_error_handler, only : fatal, warning, mom_error, mom_mesg, is_root_pe
57use mom_file_parser, only : read_param, get_param, param_file_type
60use mom_domains, only : mom_domains_init, mom_domain_type, clone_mom_domain
67
68implicit none ; private
69
72
73!>@{ CPU time clock ID
74integer :: id_clock_oda_init
75integer :: id_clock_bias_adjustment
76integer :: id_clock_apply_increments
77!>@}
78
79#include <MOM_memory.h>
80
81!> A structure with a pointer to a domain2d, to allow for the creation of arrays of pointers.
82type :: ptr_mpp_domain
83 type(domain2d), pointer :: mpp_domain => null() !< pointer to a domain2d
84end type ptr_mpp_domain
85
86!> A structure containing integer handles for bias adjustment of tracers
87type :: inc_cs
88 integer :: fldno = 0 !< The number of tracers
89 type(external_field) :: t !< The handle for the temperature file
90 type(external_field) :: s !< The handle for the salinity file
91end type inc_cs
92
93!> Control structure that contains a transpose of the ocean state across ensemble members.
94type, public :: oda_cs ; private
95 type(ocean_control_struct), pointer :: ocean_prior=> null() !< ensemble ocean prior states in DA space
96 type(ocean_control_struct), pointer :: ocean_posterior=> null() !< ensemble ocean posterior states
97 !! or increments to prior in DA space
98 type(ocean_control_struct), pointer :: ocean_increment=> null() !< A separate structure for
99 !! increment diagnostics
100 integer :: nk !< number of vertical layers used for DA
101 type(ocean_grid_type), pointer :: grid => null() !< MOM6 grid type and decomposition for the DA
102 type(ocean_grid_type), pointer :: g => null() !< MOM6 grid type and decomposition for the model
103 type(mom_domain_type), pointer, dimension(:) :: domains => null() !< Pointer to mpp_domain objects
104 !! for ensemble members
105 type(verticalgrid_type), pointer :: gv => null() !< vertical grid for DA
106 type(unit_scale_type), pointer :: &
107 us => null() !< structure containing various unit conversion factors for DA
108
109 type(domain2d), pointer :: mpp_domain => null() !< Pointer to a mpp domain object for DA
110 type(grid_type), pointer :: oda_grid !< local tracer grid
111 real, pointer, dimension(:,:,:) :: h => null() !<layer thicknesses [H ~> m or kg m-2] for DA
112 real, pointer, dimension(:,:,:) :: t_tend => null() !<layer temperature tendency from DA [C T-1 ~> degC s-1]
113 real, pointer, dimension(:,:,:) :: s_tend => null() !<layer salinity tendency from DA [S T-1 ~> ppt s-1]
114 real, pointer, dimension(:,:,:) :: t_bc_tend => null() !< The layer temperature tendency due
115 !! to bias adjustment [C T-1 ~> degC s-1]
116 real, pointer, dimension(:,:,:) :: s_bc_tend => null() !< The layer salinity tendency due
117 !! to bias adjustment [S T-1 ~> ppt s-1]
118 integer :: ni !< global i-direction grid size
119 integer :: nj !< global j-direction grid size
120 logical :: reentrant_x !< grid is reentrant in the x direction
121 logical :: reentrant_y !< grid is reentrant in the y direction
122 logical :: tripolar_n !< grid is folded at its north edge
123 logical :: symmetric !< Values at C-grid locations are symmetric
124 logical :: use_basin_mask !< If true, use a basin file to delineate weakly coupled ocean basins
125 logical :: do_bias_adjustment !< If true, use spatio-temporally varying climatological tendency
126 !! adjustment for Temperature and Salinity
127 real :: bias_adjustment_multiplier !< A scaling for the bias adjustment [nondim]
128 integer :: assim_method !< Method: NO_ASSIM,EAKF_ASSIM or OI_ASSIM
129 integer :: ensemble_size !< Size of the ensemble
130 integer :: ensemble_id = 0 !< id of the current ensemble member
131 integer, pointer, dimension(:,:) :: ensemble_pelist !< PE list for ensemble members
132 integer, pointer, dimension(:) :: filter_pelist !< PE list for ensemble members
133 real :: assim_interval !< analysis interval [T ~> s]
134 ! Profiles local to the analysis domain
135 type(ocean_profile_type), pointer :: profiles => null() !< pointer to linked list of all available profiles
136 type(ocean_profile_type), pointer :: cprofiles => null()!< pointer to linked list of current profiles
137 type(kd_root), pointer :: kdroot => null() !< A structure for storing nearest neighbors
138 type(ale_cs), pointer :: ale_cs=>null() !< ALE control structure for DA
139 logical :: use_ale_algorithm !< true is using ALE remapping
140 type(regridding_cs) :: regridcs !< ALE control structure for regridding
141 type(remapping_cs) :: remapcs !< ALE control structure for remapping
142 type(time_type) :: time !< Current Analysis time
143 type(diag_ctrl), pointer :: diag_cs=> null() !<Pointer to diagnostics control structure
144 type(inc_cs) :: inc_cs !< A Structure containing integer file handles for bias adjustment
145 integer :: id_inc_t !< A diagnostic handle for the temperature climatological adjustment
146 integer :: id_inc_s !< A diagnostic handle for the salinity climatological adjustment
147 integer :: answer_date !< The vintage of the order of arithmetic and expressions in the
148 !! remapping invoked by the ODA driver. Values below 20190101 recover
149 !! the answers from the end of 2018, while higher values use updated
150 !! and more robust forms of the same expressions.
151 logical :: reproduce_2018_nmme !< true if reproducing older NMME answers.
152end type oda_cs
153
154
155!>@{ DA parameters
156integer, parameter :: no_assim = 0, oi_assim=1, eakf_assim=2
157!>@}
158character(len=40) :: mdl = "MOM_oda_driver" !< This module's name.
159
160contains
161
162!> initialize First_guess (prior) and Analysis grid
163!! information for all ensemble members
164subroutine init_oda(Time, G, GV, US, diag_CS, CS)
165
166 type(time_type), intent(in) :: time !< The current model time.
167 type(ocean_grid_type), pointer :: g !< domain and grid information for ocean model
168 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
169 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
170 type(diag_ctrl), target, intent(inout) :: diag_cs !< A pointer to a diagnostic control structure
171 type(oda_cs), pointer, intent(inout) :: cs !< The DA control structure
172
173! Local variables
174 type(thermo_var_ptrs) :: tv_dummy
175 type(dyn_horgrid_type), pointer :: dg=> null()
176 type(hor_index_type), pointer :: hi=> null()
177 type(directories) :: dirs
178
179 type(grid_type), pointer :: t_grid !< global tracer grid
180 type(param_file_type) :: pf
181 integer :: n
182 integer :: isd, ied, jsd, jed
183 integer :: is_oda, ie_oda, js_oda, je_oda
184 integer :: isd_oda, ied_oda, jsd_oda, jed_oda
185 integer, dimension(4) :: fld_sz
186 character(len=32) :: assim_method
187 integer :: npes_pm, ens_info(6)
188 character(len=30) :: coord_mode
189 character(len=200) :: inputdir, basin_file
190 character(len=80) :: basin_var
191 character(len=80) :: remap_scheme
192 character(len=80) :: bias_correction_file, inc_file
193 integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.
194 logical :: om4_remap_via_sub_cells ! If true, use the OM4 remapping algorithm
195 real :: h_neglect, h_neglect_edge ! small thicknesses [H ~> m or kg m-2]
196
197 if (associated(cs)) call mom_error(fatal, 'Calling oda_init with associated control structure')
198 allocate(cs)
199
200 id_clock_oda_init=cpu_clock_id('(ODA initialization)')
201 call cpu_clock_begin(id_clock_oda_init)
202
203! Use ens1 parameters , this could be changed at a later time
204! if it were desirable to have alternate parameters, e.g. for the grid
205! for the analysis
206 call get_mom_input(pf,dirs,ensemble_num=0)
207 call unit_scaling_init(pf, cs%US)
208
209 call get_param(pf, mdl, "ASSIM_METHOD", assim_method, &
210 "String which determines the data assimilation method "//&
211 "Valid methods are: \'EAKF\',\'OI\', and \'NO_ASSIM\'", default='NO_ASSIM')
212 call get_param(pf, mdl, "ASSIM_INTERVAL", cs%assim_interval, &
213 "data assimilation update interval in hours",default=-1.0,units="hours",scale=3600.*us%s_to_T)
214 if (cs%assim_interval < 0.) then
215 call get_param(pf, mdl, "ASSIM_FREQUENCY", cs%assim_interval, &
216 "data assimilation update in hours. This parameter name will \n"//&
217 "be deprecated in the future. ASSIM_INTERVAL should be used instead.",default=-1.0, &
218 units="hours",scale=3600.*us%s_to_T)
219 endif
220
221 call get_param(pf, mdl, "USE_REGRIDDING", cs%use_ALE_algorithm , &
222 "If True, use the ALE algorithm (regridding/remapping).\n"//&
223 "If False, use the layered isopycnal algorithm.", default=.false. )
224 call get_param(pf, mdl, "REENTRANT_X", cs%reentrant_x, &
225 "If true, the domain is zonally reentrant.", default=.true.)
226 call get_param(pf, mdl, "REENTRANT_Y", cs%reentrant_y, &
227 "If true, the domain is meridionally reentrant.", &
228 default=.false.)
229 call get_param(pf, mdl, "TRIPOLAR_N", cs%tripolar_N, &
230 "Use tripolar connectivity at the northern edge of the "//&
231 "domain. With TRIPOLAR_N, NIGLOBAL must be even.", &
232 default=.false.)
233 call get_param(pf, mdl, "APPLY_TRACER_TENDENCY_ADJUSTMENT", cs%do_bias_adjustment, &
234 "If true, add a spatio-temporally varying climatological adjustment "//&
235 "to temperature and salinity.", &
236 default=.false.)
237 if (cs%do_bias_adjustment) then
238 call get_param(pf, mdl, "TRACER_ADJUSTMENT_FACTOR", cs%bias_adjustment_multiplier, &
239 "A multiplicative scaling factor for the climatological tracer tendency adjustment ", &
240 units="nondim", default=1.0)
241 endif
242 call get_param(pf, mdl, "USE_BASIN_MASK", cs%use_basin_mask, &
243 "If true, add a basin mask to delineate weakly connected "//&
244 "ocean basins for the purpose of data assimilation.", &
245 default=.false.)
246
247 call get_param(pf, mdl, "NIGLOBAL", cs%ni, &
248 "The total number of thickness grid points in the "//&
249 "x-direction in the physical domain.")
250 call get_param(pf, mdl, "NJGLOBAL", cs%nj, &
251 "The total number of thickness grid points in the "//&
252 "y-direction in the physical domain.")
253 call get_param(pf, mdl, "INPUTDIR", inputdir)
254 call get_param(pf, mdl, "ODA_REMAPPING_SCHEME", remap_scheme, &
255 "This sets the reconstruction scheme used "//&
256 "for vertical remapping for all ODA variables. "//&
257 "It can be one of the following schemes: "//&
258 trim(remappingschemesdoc), default="PPM_H4")
259 call get_param(pf, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
260 "This sets the default value for the various _ANSWER_DATE parameters.", &
261 default=99991231)
262 call get_param(pf, mdl, "ODA_ANSWER_DATE", cs%answer_date, &
263 "The vintage of the order of arithmetic and expressions used by the ODA driver "//&
264 "Values below 20190101 recover the answers from the end of 2018, while higher "//&
265 "values use updated and more robust forms of the same expressions.", &
266 default=default_answer_date, do_not_log=.not.gv%Boussinesq)
267 if (.not.gv%Boussinesq) cs%answer_date = max(cs%answer_date, 20230701)
268
269 call get_param(pf, mdl, "REPRODUCE_2018_NMME_ANSWERS", cs%reproduce_2018_nmme, &
270 "Logical flag needed to reproduce older NMME forecast answers. "//&
271 "True gives old answers, the default of false gives different answers.", &
272 default=.false.)
273
274 inputdir = slasher(inputdir)
275
276 select case(lowercase(trim(assim_method)))
277 case('eakf')
278 cs%assim_method = eakf_assim
279 case('oi')
280 cs%assim_method = oi_assim
281 case('no_assim')
282 cs%assim_method = no_assim
283 case default
284 call mom_error(fatal, "Invalid assimilation method provided")
285 end select
286
287 ens_info = get_ensemble_size()
288 cs%ensemble_size = ens_info(1)
289 npes_pm=ens_info(3)
290 cs%ensemble_id = get_ensemble_id()
291 !! Switch to global pelist
292 allocate(cs%ensemble_pelist(cs%ensemble_size,npes_pm))
293 allocate(cs%filter_pelist(cs%ensemble_size*npes_pm))
294 call get_ensemble_pelist(cs%ensemble_pelist, 'ocean')
295 call get_ensemble_filter_pelist(cs%filter_pelist, 'ocean')
296
297 call set_pelist(cs%filter_pelist)
298
299 allocate(cs%domains(cs%ensemble_size))
300 cs%domains(cs%ensemble_id)%mpp_domain => g%Domain%mpp_domain ! this should go away
301 do n=1,cs%ensemble_size
302 if (.not. associated(cs%domains(n)%mpp_domain)) allocate(cs%domains(n)%mpp_domain)
303 call set_rootpe(cs%ensemble_pelist(n,1)) ! this line is not in Feiyu's version (needed?)
304 call broadcast_domain(cs%domains(n)%mpp_domain)
305 enddo
306 call set_rootpe(cs%filter_pelist(1)) ! this line is not in Feiyu's version (needed?)
307 cs%G => g
308 allocate(cs%Grid)
309 ! params NIHALO_ODA, NJHALO_ODA set the DA halo size
310 call mom_domains_init(cs%Grid%Domain, pf, param_suffix='_ODA', us=cs%US)
311 allocate(hi)
312 call hor_index_init(cs%Grid%Domain, hi, pf)
313 call verticalgridinit( pf, cs%GV, cs%US )
314 allocate(dg)
315 call create_dyn_horgrid(dg, hi)
316 call clone_mom_domain(cs%Grid%Domain, dg%Domain,symmetric=.false.)
317 call set_grid_metrics(dg, pf, cs%US)
318 call mom_initialize_topography(dg%bathyT, dg%max_depth, dg, pf, cs%US)
319 call mom_initialize_coord(cs%GV, cs%US, pf, tv_dummy, dg%max_depth)
320 call ale_init(pf, cs%G, cs%GV, cs%US, dg%max_depth, cs%ALE_CS)
321 call mom_grid_init(cs%Grid, pf, global_indexing=.false.)
322 call ale_updateverticalgridtype(cs%ALE_CS, cs%GV)
323 call copy_dyngrid_to_mom_grid(dg, cs%Grid, cs%US)
324 cs%mpp_domain => cs%Grid%Domain%mpp_domain
325 cs%Grid%ke = cs%GV%ke
326 cs%nk = cs%GV%ke
327 ! initialize storage for prior and posterior
328 allocate(cs%Ocean_prior)
329 call init_ocean_ensemble(cs%Ocean_prior,cs%Grid,cs%GV,cs%ensemble_size)
330 allocate(cs%Ocean_posterior)
331 call init_ocean_ensemble(cs%Ocean_posterior,cs%Grid,cs%GV,cs%ensemble_size)
332 allocate(cs%Ocean_increment)
333 call init_ocean_ensemble(cs%Ocean_increment,cs%Grid,cs%GV,cs%ensemble_size)
334
335
336 call get_param(pf, 'oda_driver', "REGRIDDING_COORDINATE_MODE", coord_mode, &
337 "Coordinate mode for vertical regridding.", &
338 default="ZSTAR", fail_if_missing=.false.)
339 call get_param(pf, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
340 do_not_log=.true., default=.true.)
341 call get_param(pf, mdl, "ODA_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
342 "If true, use the OM4 remapping-via-subcells algorithm for ODA. "//&
343 "See REMAPPING_USE_OM4_SUBCELLS for more details. "//&
344 "We recommend setting this option to false.", default=om4_remap_via_sub_cells)
345 call initialize_regridding(cs%regridCS, cs%G, cs%GV, cs%US, dg%max_depth,pf,'oda_driver',coord_mode,'','')
346
347 h_neglect = set_h_neglect(gv, cs%answer_date, h_neglect_edge)
348 call initialize_remapping(cs%remapCS, remap_scheme, om4_remap_via_sub_cells=om4_remap_via_sub_cells, &
349 h_neglect=h_neglect, h_neglect_edge=h_neglect_edge, answer_date=cs%answer_date)
350 call set_regrid_params(cs%regridCS, min_thickness=0.)
351 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
352
353 ! breaking with the MOM6 convention and using global indices
354 !call get_domain_extent(G%Domain,is,ie,js,je,isd,ied,jsd,jed,&
355 ! isg,ieg,jsg,jeg,idg_offset,jdg_offset,symmetric)
356 !isd = isd+idg_offset ; ied = ied+idg_offset ! using global indexing within the DA module
357 !jsd = jsd+jdg_offset ; jed = jed+jdg_offset ! TODO: switch to local indexing? (mjh)
358
359 if (.not. associated(cs%h)) then
360 allocate(cs%h(isd:ied,jsd:jed,cs%GV%ke), source=cs%GV%Angstrom_H)
361 ! assign thicknesses
362 call ale_initthicknesstocoord(cs%ALE_CS, g, cs%GV, cs%h)
363 endif
364
365 allocate(cs%T_tend(isd:ied,jsd:jed,cs%GV%ke), source=0.0)
366 allocate(cs%S_tend(isd:ied,jsd:jed,cs%GV%ke), source=0.0)
367! call set_axes_info(CS%Grid, CS%GV, CS%US, PF, CS%diag_cs, set_vertical=.true.) ! missing in Feiyu's fork
368 allocate(cs%oda_grid)
369 cs%oda_grid%x => cs%Grid%geolonT
370 cs%oda_grid%y => cs%Grid%geolatT
371
372
373 if (cs%use_basin_mask) then
374 call get_param(pf, 'oda_driver', "BASIN_FILE", basin_file, &
375 "A file in which to find the basin masks.", default="basin.nc")
376 basin_file = trim(inputdir) // trim(basin_file)
377 call get_param(pf, 'oda_driver', "BASIN_VAR", basin_var, &
378 "The basin mask variable in BASIN_FILE.", default="basin")
379 ! Need different data domain indices for the ODA ensemble basin mask.
380 call get_domain_extent(cs%Grid%Domain, is_oda, ie_oda, js_oda, je_oda, isd_oda, ied_oda, jsd_oda, jed_oda)
381 allocate(cs%oda_grid%basin_mask(isd_oda:ied_oda,jsd_oda:jed_oda), source=0.0)
382 call mom_read_data(basin_file, basin_var, cs%oda_grid%basin_mask, cs%Grid%domain, timelevel=1)
383 endif
384
385 ! set up diag variables for analysis increments
386 cs%diag_CS => diag_cs
387 cs%id_inc_t = register_diag_field('ocean_model', 'temp_increment', diag_cs%axesTL, &
388 time, 'ocean potential temperature increments', 'degC', conversion=us%C_to_degC)
389 cs%id_inc_s = register_diag_field('ocean_model', 'salt_increment', diag_cs%axesTL, &
390 time, 'ocean salinity increments', 'psu', conversion=us%S_to_ppt)
391
392 !! get global grid information from ocean model needed for ODA initialization
393 t_grid => null()
394 call set_up_global_tgrid(t_grid, cs, g)
395
396 call ocean_da_core_init(cs%mpp_domain, t_grid, cs%Profiles, time)
397 deallocate(t_grid)
398 cs%Time = time
399 !! switch back to ensemble member pelist
400 call set_pelist(cs%ensemble_pelist(cs%ensemble_id,:))
401
402 if (cs%do_bias_adjustment) then
403 call get_param(pf, mdl, "TEMP_SALT_ADJUSTMENT_FILE", bias_correction_file, &
404 "The name of the file containing temperature and salinity "//&
405 "tendency adjustments", default='temp_salt_adjustment.nc')
406
407 inc_file = trim(inputdir) // trim(bias_correction_file)
408 cs%INC_CS%T = init_extern_field(inc_file, "temp_increment", &
409 correct_leap_year_inconsistency=.true.,verbose=.true.,domain=g%Domain%mpp_domain)
410 cs%INC_CS%S = init_extern_field(inc_file, "salt_increment", &
411 correct_leap_year_inconsistency=.true.,verbose=.true.,domain=g%Domain%mpp_domain)
412 call get_external_field_info(cs%INC_CS%T, size=fld_sz)
413 cs%INC_CS%fldno = 2
414 if (cs%nk /= fld_sz(3)) call mom_error(fatal,'Increment levels /= ODA levels')
415
416 allocate(cs%T_bc_tend(g%isd:g%ied,g%jsd:g%jed,cs%GV%ke), source=0.0)
417 allocate(cs%S_bc_tend(g%isd:g%ied,g%jsd:g%jed,cs%GV%ke), source=0.0)
418 endif
419
420 call cpu_clock_end(id_clock_oda_init)
421
422! if (CS%write_obs) then
423! temp_fid = open_profile_file("temp_"//trim(obs_file))
424! salt_fid = open_profile_file("salt_"//trim(obs_file))
425! endif
426
427end subroutine init_oda
428
429!> Copy ensemble member tracers to ensemble vector.
430subroutine set_prior_tracer(Time, G, GV, h, tv, CS)
431 type(time_type), intent(in) :: time !< The current model time
432 type(ocean_grid_type), pointer :: g !< domain and grid information for ocean model
433 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
434 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
435 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables
436
437 type(oda_cs), pointer :: cs !< ocean DA control structure
438 real, dimension(SZI_(G),SZJ_(G),CS%nk) :: t ! Temperature on the analysis grid [C ~> degC]
439 real, dimension(SZI_(G),SZJ_(G),CS%nk) :: s ! Salinity on the analysis grid [S ~> ppt]
440 integer :: i, j, m
441 integer :: isc, iec, jsc, jec
442
443 ! return if not time for analysis
444 if (time < cs%Time) return
445
446 if (.not. associated(cs%Grid)) call mom_error(fatal,'ODA_CS ensemble horizontal grid not associated')
447 if (.not. associated(cs%GV)) call mom_error(fatal,'ODA_CS ensemble vertical grid not associated')
448
449 !! switch to global pelist
450 call set_pelist(cs%filter_pelist)
451 !call MOM_mesg('Setting prior')
452
453 ! computational domain for the analysis grid
454 isc = cs%Grid%isc ; iec = cs%Grid%iec ; jsc = cs%Grid%jsc ; jec = cs%Grid%jec
455 ! array extents for the ensemble member
456 !call get_domain_extent(CS%domains(CS%ensemble_id),is,ie,js,je,isd,ied,jsd,jed,&
457 ! isg,ieg,jsg,jeg,idg_offset,jdg_offset,symmetric)
458 ! remap temperature and salinity from the ensemble member to the analysis grid
459 do j=g%jsc,g%jec ; do i=g%isc,g%iec
460 call remapping_core_h(cs%remapCS, gv%ke, h(i,j,:), tv%T(i,j,:), &
461 cs%nk, cs%h(i,j,:), t(i,j,:))
462 call remapping_core_h(cs%remapCS, gv%ke, h(i,j,:), tv%S(i,j,:), &
463 cs%nk, cs%h(i,j,:), s(i,j,:))
464 enddo ; enddo
465 ! cast ensemble members to the analysis domain
466 do m=1,cs%ensemble_size
467 call redistribute_array(cs%domains(m)%mpp_domain, t,&
468 cs%mpp_domain, cs%Ocean_prior%T(:,:,:,m), complete=.true.)
469 call redistribute_array(cs%domains(m)%mpp_domain, s,&
470 cs%mpp_domain, cs%Ocean_prior%S(:,:,:,m), complete=.true.)
471 enddo
472
473 do m=1,cs%ensemble_size
474 call pass_var(cs%Ocean_prior%T(:,:,:,m),cs%Grid%domain)
475 call pass_var(cs%Ocean_prior%S(:,:,:,m),cs%Grid%domain)
476 enddo
477
478 !! switch back to ensemble member pelist
479 call set_pelist(cs%ensemble_pelist(cs%ensemble_id,:))
480
481 return
482
483end subroutine set_prior_tracer
484
485!> Returns posterior adjustments or full state
486!!Note that only those PEs associated with an ensemble member receive data
487subroutine get_posterior_tracer(Time, CS, increment)
488 type(time_type), intent(in) :: time !< the current model time
489 type(oda_cs), pointer :: cs !< ocean DA control structure
490 logical, optional, intent(in) :: increment !< True if returning increment only
491
492 type(ocean_control_struct), pointer :: ocean_increment=>null()
493 integer :: m
494 logical :: get_inc
495
496
497 ! return if not analysis time (retain pointers for h and tv)
498 if (time < cs%Time .or. cs%assim_method == no_assim) return
499
500
501 !! switch to global pelist
502 call set_pelist(cs%filter_pelist)
503 call mom_mesg('Getting posterior')
504
505 !! Calculate and redistribute increments to CS%tv right after assimilation
506 !! Retain CS%tv to calculate increments for IAU updates CS%tv_inc otherwise
507 get_inc = .true.
508 if (present(increment)) get_inc = increment
509
510 if (get_inc) then
511 cs%Ocean_increment%T = cs%Ocean_posterior%T - cs%Ocean_prior%T
512 cs%Ocean_increment%S = cs%Ocean_posterior%S - cs%Ocean_prior%S
513 endif
514 ! It may be necessary to check whether the increment and ocean state have the
515 ! same dimensionally rescaled units.
516 do m=1,cs%ensemble_size
517 if (get_inc) then
518 call redistribute_array(cs%mpp_domain, cs%Ocean_increment%T(:,:,:,m),&
519 cs%domains(m)%mpp_domain, cs%T_tend, complete=.true.)
520 call redistribute_array(cs%mpp_domain, cs%Ocean_increment%S(:,:,:,m),&
521 cs%domains(m)%mpp_domain, cs%S_tend, complete=.true.)
522 else
523 call redistribute_array(cs%mpp_domain, cs%Ocean_posterior%T(:,:,:,m),&
524 cs%domains(m)%mpp_domain, cs%T_tend, complete=.true.)
525 call redistribute_array(cs%mpp_domain, cs%Ocean_posterior%S(:,:,:,m),&
526 cs%domains(m)%mpp_domain, cs%S_tend, complete=.true.)
527 endif
528 enddo
529
530
531 !! switch back to ensemble member pelist
532 call set_pelist(cs%ensemble_pelist(cs%ensemble_id,:))
533
534 call pass_var(cs%T_tend,cs%domains(cs%ensemble_id))
535 call pass_var(cs%S_tend,cs%domains(cs%ensemble_id))
536
537 !convert to a tendency (degC or PSU per second)
538 cs%T_tend = cs%T_tend / (cs%assim_interval)
539 cs%S_tend = cs%S_tend / (cs%assim_interval)
540
541
542end subroutine get_posterior_tracer
543
544!> Gather observations and call ODA routines
545subroutine oda(Time, CS)
546 type(time_type), intent(in) :: time !< the current model time
547 type(oda_cs), pointer :: cs !< A pointer the ocean DA control structure
548
549 if ( time >= cs%Time ) then
550
551 !! switch to global pelist
552 call set_pelist(cs%filter_pelist)
553 call get_profiles(time, cs%Profiles, cs%CProfiles)
554#ifdef ENABLE_ECDA
555 call ensemble_filter(cs%Ocean_prior, cs%Ocean_posterior, cs%CProfiles, cs%kdroot, cs%mpp_domain, cs%oda_grid)
556#endif
557 !! switch back to ensemble member pelist
558 call set_pelist(cs%ensemble_pelist(cs%ensemble_id,:))
559 call get_posterior_tracer(time, cs, increment=.true.)
560 if (cs%do_bias_adjustment) call get_bias_correction_tracer(time, cs%US, cs)
561
562 endif
563
564 return
565end subroutine oda
566
567subroutine get_bias_correction_tracer(Time, US, CS)
568 type(time_type), intent(in) :: Time !< the current model time
569 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
570 type(oda_cs), pointer :: CS !< ocean DA control structure
571
572 ! Local variables
573 real, allocatable, dimension(:,:,:) :: T_bias ! Estimated temperature tendency bias [C T-1 ~> degC s-1]
574 real, allocatable, dimension(:,:,:) :: S_bias ! Estimated salinity tendency bias [S T-1 ~> ppt s-1]
575 real, allocatable, dimension(:,:,:) :: valid_flag ! Valid value flag on the horizontal model grid
576 ! and input-file vertical levels [nondim]
577 real, allocatable, dimension(:), target :: z_in ! Cell center depths for input data [Z ~> m]
578 real, allocatable, dimension(:), target :: z_edges_in ! Cell edge depths for input data [Z ~> m]
579 real :: missing_value ! A value indicating that there is no valid input data at this point [CU ~> conc]
580 integer, dimension(3) :: fld_sz
581 integer :: i,j,k
582
583
584 call cpu_clock_begin(id_clock_bias_adjustment)
585 call horiz_interp_and_extrap_tracer(cs%INC_CS%T, time, cs%G, t_bias, &
586 valid_flag, z_in, z_edges_in, missing_value, scale=us%degC_to_C*us%s_to_T, spongeongrid=.true., &
587 answer_date=cs%answer_date)
588 call horiz_interp_and_extrap_tracer(cs%INC_CS%S, time, cs%G, s_bias, &
589 valid_flag, z_in, z_edges_in, missing_value, scale=us%ppt_to_S*us%s_to_T, spongeongrid=.true., &
590 answer_date=cs%answer_date)
591
592 ! This should be replaced to use mask_z instead of the following lines
593 ! which are intended to zero land values using an arbitrary limit.
594 fld_sz=shape(t_bias)
595 if (cs%reproduce_2018_nmme) then
596 do i=1,fld_sz(1)
597 do j=1,fld_sz(2)
598 do k=1,fld_sz(3)
599 ! The following two lines are needed for backward compatibility for NMME answers (2018 vintage)
600 ! These were implemented to catch missing values, so large values are excluded.
601 if (t_bias(i,j,k) > 1.0e-3*us%degC_to_C) t_bias(i,j,k) = 0.0
602 if (s_bias(i,j,k) > 1.0e-3*us%ppt_to_S) s_bias(i,j,k) = 0.0
603 enddo
604 enddo
605 enddo
606 else
607 do i=1,fld_sz(1)
608 do j=1,fld_sz(2)
609 do k=1,fld_sz(3)
610 if (valid_flag(i,j,k)==0.) then
611 t_bias(i,j,k)=0.0
612 s_bias(i,j,k)=0.0
613 endif
614 enddo
615 enddo
616 enddo
617 endif
618
619 cs%T_bc_tend = t_bias * cs%bias_adjustment_multiplier
620 cs%S_bc_tend = s_bias * cs%bias_adjustment_multiplier
621
622 call pass_var(cs%T_bc_tend, cs%domains(cs%ensemble_id))
623 call pass_var(cs%S_bc_tend, cs%domains(cs%ensemble_id))
624
625 call cpu_clock_end(id_clock_bias_adjustment)
626
627end subroutine get_bias_correction_tracer
628
629!> Finalize DA module
630subroutine oda_end(CS)
631 type(oda_cs), intent(inout) :: cs !< the ocean DA control structure
632
633end subroutine oda_end
634
635!> Initialize DA module
636subroutine init_ocean_ensemble(CS, Grid, GV, ens_size)
637 type(ocean_control_struct), pointer :: CS !< Pointer to ODA control structure
638 type(ocean_grid_type), pointer :: Grid !< Pointer to ocean analysis grid
639 type(verticalgrid_type), pointer :: GV !< Pointer to DA vertical grid
640 integer, intent(in) :: ens_size !< ensemble size
641
642 integer :: is, ie, js, je, nk
643
644 nk = gv%ke
645 is = grid%isd ; ie = grid%ied
646 js = grid%jsd ; je = grid%jed
647 cs%ensemble_size = ens_size
648 allocate(cs%T(is:ie,js:je,nk,ens_size))
649 allocate(cs%S(is:ie,js:je,nk,ens_size))
650 allocate(cs%SSH(is:ie,js:je,ens_size))
651! allocate(CS%id_t(ens_size), source=-1)
652! allocate(CS%id_s(ens_size), source=-1)
653! allocate(CS%U(is:ie,js:je,nk,ens_size))
654! allocate(CS%V(is:ie,js:je,nk,ens_size))
655! allocate(CS%id_u(ens_size), source=-1)
656! allocate(CS%id_v(ens_size), source=-1)
657! allocate(CS%id_ssh(ens_size), source=-1)
658
659 return
660end subroutine init_ocean_ensemble
661
662!> Set the next analysis time
663subroutine set_analysis_time(Time, CS)
664 type(time_type), intent(in) :: time !< the current model time
665 type(oda_cs), pointer, intent(inout) :: cs !< the DA control structure
666
667 character(len=160) :: mesg ! The text of an error message
668 integer :: yr, mon, day, hr, min, sec
669
670 if (time >= cs%Time) then
671 ! increment the analysis time to the next step
672 cs%Time = cs%Time + real_to_time(cs%assim_interval, unscale=cs%US%T_to_s)
673
674 call get_date(time, yr, mon, day, hr, min, sec)
675 write(mesg,*) 'Model Time: ', yr, mon, day, hr, min, sec
676 call mom_mesg("set_analysis_time: "//trim(mesg))
677 call get_date(cs%time, yr, mon, day, hr, min, sec)
678 write(mesg,*) 'Assimilation Time: ', yr, mon, day, hr, min, sec
679 call mom_mesg("set_analysis_time: "//trim(mesg))
680 endif
681 if (cs%Time < time) then
682 call mom_error(fatal, " set_analysis_time: " // &
683 "assimilation interval appears to be shorter than " // &
684 "the model timestep")
685 endif
686 return
687
688end subroutine set_analysis_time
689
690
691!> Apply increments to tracers
692subroutine apply_oda_tracer_increments(dt, Time_end, G, GV, tv, h, CS)
693 real, intent(in) :: dt !< The tracer timestep [T ~> s]
694 type(time_type), intent(in) :: time_end !< Time at the end of the interval
695 type(ocean_grid_type), intent(in) :: g !< ocean grid structure
696 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
697 type(thermo_var_ptrs), intent(inout) :: tv !< A structure pointing to various thermodynamic variables
698 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
699 intent(in) :: h !< layer thickness [H ~> m or kg m-2]
700 type(oda_cs), pointer :: cs !< the data assimilation structure
701
702 !! local variables
703 integer :: i, j
704 integer :: isc, iec, jsc, jec
705 real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: t_tend_inc !< an adjustment to the temperature
706 !! tendency [C T-1 ~> degC s-1]
707 real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: s_tend_inc !< an adjustment to the salinity
708 !! tendency [S T-1 ~> ppt s-1]
709 real, dimension(SZI_(G),SZJ_(G),SZK_(CS%Grid)) :: t_tend !< The temperature tendency adjustment from
710 !! DA [C T-1 ~> degC s-1]
711 real, dimension(SZI_(G),SZJ_(G),SZK_(CS%Grid)) :: s_tend !< The salinity tendency adjustment from DA
712 !! [S T-1 ~> ppt s-1]
713
714 if (.not. associated(cs)) return
715 if (cs%assim_method == no_assim .and. (.not. cs%do_bias_adjustment)) return
716
717 call cpu_clock_begin(id_clock_apply_increments)
718
719 t_tend_inc(:,:,:) = 0.0 ; s_tend_inc(:,:,:) = 0.0 ; t_tend(:,:,:) = 0.0 ; s_tend(:,:,:) = 0.0
720 if (cs%assim_method > 0 ) then
721 t_tend = t_tend + cs%T_tend
722 s_tend = s_tend + cs%S_tend
723 endif
724 if (cs%do_bias_adjustment ) then
725 t_tend = t_tend + cs%T_bc_tend
726 s_tend = s_tend + cs%S_bc_tend
727 endif
728
729 isc=g%isc ; iec=g%iec ; jsc=g%jsc ; jec=g%jec
730 do j=jsc,jec ; do i=isc,iec
731 call remapping_core_h(cs%remapCS, cs%nk, cs%h(i,j,:), t_tend(i,j,:), &
732 g%ke, h(i,j,:), t_tend_inc(i,j,:))
733 call remapping_core_h(cs%remapCS, cs%nk, cs%h(i,j,:), s_tend(i,j,:), &
734 g%ke, h(i,j,:), s_tend_inc(i,j,:))
735 enddo ; enddo
736
737
738 call pass_var(t_tend_inc, g%Domain)
739 call pass_var(s_tend_inc, g%Domain)
740
741 tv%T(isc:iec,jsc:jec,:) = tv%T(isc:iec,jsc:jec,:) + t_tend_inc(isc:iec,jsc:jec,:)*dt
742 tv%S(isc:iec,jsc:jec,:) = tv%S(isc:iec,jsc:jec,:) + s_tend_inc(isc:iec,jsc:jec,:)*dt
743
744 call pass_var(tv%T, g%Domain)
745 call pass_var(tv%S, g%Domain)
746
747 call enable_averaging(dt, time_end, cs%diag_CS)
748 if (cs%id_inc_t > 0) call post_data(cs%id_inc_t, t_tend_inc, cs%diag_CS)
749 if (cs%id_inc_s > 0) call post_data(cs%id_inc_s, s_tend_inc, cs%diag_CS)
750 call disable_averaging(cs%diag_CS)
751
752 call diag_update_remap_grids(cs%diag_CS)
753 call cpu_clock_end(id_clock_apply_increments)
754
755
756end subroutine apply_oda_tracer_increments
757
758!> Set up the grid of thicknesses at tracer points throughout the global domain
759 subroutine set_up_global_tgrid(T_grid, CS, G)
760 type(grid_type), pointer :: T_grid !< global tracer grid
761 type(oda_cs), pointer, intent(in) :: CS !< A pointer to DA control structure.
762 type(ocean_grid_type), pointer :: G !< domain and grid information for ocean model
763
764 ! local variables
765 real, dimension(:,:), allocatable :: &
766 global2D, & ! A layer thickness in the entire global domain [H ~> m or kg m-2]
767 global2D_old ! The thickness of the layer above the one in global2D in the entire
768 ! global domain [H ~> m or kg m-2]
769 integer :: i, j, k
770
771 ! get global grid information from ocean_model
772 t_grid=>null()
773 !if (associated(T_grid)) call MOM_error(FATAL,'MOM_oda_driver:set_up_global_tgrid called with associated T_grid')
774
775 allocate(t_grid)
776 t_grid%ni = cs%ni
777 t_grid%nj = cs%nj
778 t_grid%nk = cs%nk
779 allocate(t_grid%x(cs%ni,cs%nj))
780 allocate(t_grid%y(cs%ni,cs%nj))
781 allocate(t_grid%bathyT(cs%ni,cs%nj))
782 call global_field(cs%mpp_domain, cs%Grid%geolonT, t_grid%x)
783 call global_field(cs%mpp_domain, cs%Grid%geolatT, t_grid%y)
784 call global_field(cs%domains(cs%ensemble_id)%mpp_domain, g%bathyT, t_grid%bathyT)
785 if (cs%use_basin_mask) then
786 allocate(t_grid%basin_mask(cs%ni,cs%nj))
787 call global_field(cs%mpp_domain, cs%oda_grid%basin_mask, t_grid%basin_mask)
788 endif
789 allocate(t_grid%mask(cs%ni,cs%nj,cs%nk), source=0.0)
790 allocate(t_grid%z(cs%ni,cs%nj,cs%nk), source=0.0)
791 allocate(global2d(cs%ni,cs%nj))
792 allocate(global2d_old(cs%ni,cs%nj))
793
794 do k = 1, cs%nk
795 call global_field(g%Domain%mpp_domain, cs%h(:,:,k), global2d)
796 do i=1,cs%ni ; do j=1,cs%nj
797 ! ###Does the next line need to be revised? Perhaps it should be
798 ! if ( global2D(i,j) > 1.0*GV%H_to_m ) then
799 if ( global2d(i,j) > 1 ) then
800 t_grid%mask(i,j,k) = 1.0
801 endif
802 enddo ; enddo
803 if (k == 1) then
804 t_grid%z(:,:,k) = global2d/2
805 else
806 t_grid%z(:,:,k) = t_grid%z(:,:,k-1) + (global2d + global2d_old)/2
807 endif
808 global2d_old = global2d
809 enddo
810
811 deallocate(global2d)
812 deallocate(global2d_old)
813 end subroutine set_up_global_tgrid
814
815!> \namespace MOM_oda_driver_mod
816!!
817!! \section section_ODA The Ocean data assimilation (DA) and Ensemble Framework
818!!
819!! The DA framework implements ensemble capability in MOM6. Currently, this framework
820!! is enabled using the cpp directive ENSEMBLE_OCEAN. The ensembles need to be generated
821!! at the level of the calling routine for oda_init or above. The ensemble instances may
822!! exist on overlapping or non-overlapping processors. The ensemble information is accessed
823!! via the FMS ensemble manager. An independent PE layout is used to gather (prior) ensemble
824!! member information where this information is stored in the ODA control structure. This
825!! module was developed in collaboration with Feiyu Lu and Tony Rosati in the GFDL prediction
826!! group for use in their coupled ensemble framework. These interfaces should be suitable for
827!! interfacing MOM6 to other data assimilation packages as well.
828
829end module mom_oda_driver_mod