MOM6
MOM_oda_driver.F90
Go to the documentation of this file.
1 !> Interfaces for MOM6 ensembles and data assimilation.
3 
4 ! This file is part of MOM6. see LICENSE.md for the license.
5 use fms_mod, only : open_namelist_file, close_file, check_nml_error
6 use fms_mod, only : error_mesg, fatal
7 use mpp_mod, only : stdout, stdlog, mpp_error, npes=>mpp_npes,pe=>mpp_pe
8 use mpp_mod, only : set_current_pelist => mpp_set_current_pelist
9 use mpp_mod, only : set_root_pe => mpp_set_root_pe
10 use mpp_mod, only : mpp_sync_self, mpp_sum, get_pelist=>mpp_get_current_pelist, mpp_root_pe
11 use mpp_mod, only : set_stack_size=>mpp_set_stack_size, broadcast=>mpp_broadcast
12 use mpp_io_mod, only : io_set_stack_size=>mpp_io_set_stack_size
13 use mpp_io_mod, only : mpp_single,mpp_multi
14 use mpp_domains_mod, only : domain2d, mpp_global_field
15 use mpp_domains_mod, only : mpp_get_compute_domain, mpp_get_data_domain
16 use mpp_domains_mod, only : mpp_redistribute, mpp_broadcast_domain
17 use mpp_domains_mod, only : set_domains_stack_size=>mpp_domains_set_stack_size
18 use diag_manager_mod, only : register_diag_field, diag_axis_init, send_data
19 use ensemble_manager_mod, only : get_ensemble_id, get_ensemble_size
20 use ensemble_manager_mod, only : get_ensemble_pelist, get_ensemble_filter_pelist
21 use time_manager_mod, only : time_type, decrement_time, increment_time
22 use time_manager_mod, only : get_date, operator(>=),operator(/=),operator(==),operator(<)
23 use constants_mod, only : radius, epsln
24 ! ODA Modules
25 use ocean_da_types_mod, only : grid_type, ocean_profile_type, ocean_control_struct
26 use ocean_da_core_mod, only : ocean_da_core_init, get_profiles
27 !use eakf_oda_mod, only : ensemble_filter
28 use write_ocean_obs_mod, only : open_profile_file
29 use write_ocean_obs_mod, only : write_profile,close_profile_file
30 use kdtree, only : kd_root !# JEDI
31 ! MOM Modules
32 use mom_io, only : slasher, mom_read_data
34 use mom_error_handler, only : fatal, warning, mom_error, mom_mesg, is_root_pe
53 
54 implicit none ; private
55 
58 
59 #include <MOM_memory.h>
60 
61 !> Control structure that contains a transpose of the ocean state across ensemble members.
62 type, public :: oda_cs ; private
63  type(ocean_control_struct), pointer :: ocean_prior=> null() !< ensemble ocean prior states in DA space
64  type(ocean_control_struct), pointer :: ocean_posterior=> null() !< ensemble ocean posterior states
65  !! or increments to prior in DA space
66  integer :: nk !< number of vertical layers used for DA
67  type(ocean_grid_type), pointer :: grid => null() !< MOM6 grid type and decomposition for the DA
68  type(ptr_mpp_domain), pointer, dimension(:) :: domains => null() !< Pointer to mpp_domain objects
69  !! for ensemble members
70  type(verticalgrid_type), pointer :: gv => null() !< vertical grid for DA
71  type(unit_scale_type), pointer :: &
72  us => null() !< structure containing various unit conversion factors for DA
73 
74  type(domain2d), pointer :: mpp_domain => null() !< Pointer to a mpp domain object for DA
75  type(grid_type), pointer :: oda_grid !< local tracer grid
76  real, pointer, dimension(:,:,:) :: h => null() !<layer thicknesses [H ~> m or kg m-2] for DA
77  type(thermo_var_ptrs), pointer :: tv => null() !< pointer to thermodynamic variables
78  integer :: ni !< global i-direction grid size
79  integer :: nj !< global j-direction grid size
80  logical :: reentrant_x !< grid is reentrant in the x direction
81  logical :: reentrant_y !< grid is reentrant in the y direction
82  logical :: tripolar_n !< grid is folded at its north edge
83  logical :: symmetric !< Values at C-grid locations are symmetric
84  integer :: assim_method !< Method: NO_ASSIM,EAKF_ASSIM or OI_ASSIM
85  integer :: ensemble_size !< Size of the ensemble
86  integer :: ensemble_id = 0 !< id of the current ensemble member
87  integer, pointer, dimension(:,:) :: ensemble_pelist !< PE list for ensemble members
88  integer, pointer, dimension(:) :: filter_pelist !< PE list for ensemble members
89  integer :: assim_frequency !< analysis interval in hours
90  ! Profiles local to the analysis domain
91  type(ocean_profile_type), pointer :: profiles => null() !< pointer to linked list of all available profiles
92  type(ocean_profile_type), pointer :: cprofiles => null()!< pointer to linked list of current profiles
93  type(kd_root), pointer :: kdroot => null() !< A structure for storing nearest neighbors
94  type(ale_cs), pointer :: ale_cs=>null() !< ALE control structure for DA
95  logical :: use_ale_algorithm !< true is using ALE remapping
96  type(regridding_cs) :: regridcs !< ALE control structure for regridding
97  type(remapping_cs) :: remapcs !< ALE control structure for remapping
98  type(time_type) :: time !< Current Analysis time
99  type(diag_ctrl) :: diag_cs !<Diagnostics control structure
100 end type oda_cs
101 
102 !> A structure with a pointer to a domain2d, to allow for the creation of arrays of pointers.
104  type(domain2d), pointer :: mpp_domain => null() !< pointer to an mpp domain2d
105 end type ptr_mpp_domain
106 
107 !>@{ DA parameters
108 integer, parameter :: no_assim = 0, oi_assim=1, eakf_assim=2
109 !!@}
110 
111 contains
112 
113 !> initialize First_guess (prior) and Analysis grid
114 !! information for all ensemble members
115 subroutine init_oda(Time, G, GV, CS)
116 
117  type(time_type), intent(in) :: time !< The current model time.
118  type(ocean_grid_type), pointer :: g !< domain and grid information for ocean model
119  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
120  type(oda_cs), pointer, intent(inout) :: cs !< The DA control structure
121 
122 ! Local variables
123  type(thermo_var_ptrs) :: tv_dummy
124  type(dyn_horgrid_type), pointer :: dg=> null()
125  type(hor_index_type), pointer :: hi=> null()
126  type(directories) :: dirs
127 
128  type(grid_type), pointer :: t_grid !< global tracer grid
129  real, dimension(:,:), allocatable :: global2d, global2d_old
130  real, dimension(:), allocatable :: lon1d, lat1d, glon1d, glat1d
131  type(param_file_type) :: pf
132  integer :: n, m, k, i, j, nk
133  integer :: is,ie,js,je,isd,ied,jsd,jed
134  integer :: stdout_unit
135  character(len=32) :: assim_method
136  integer :: npes_pm, ens_info(6), ni, nj
137  character(len=128) :: mesg
138  character(len=32) :: fldnam
139  character(len=30) :: coord_mode
140  character(len=200) :: inputdir, basin_file
141  logical :: reentrant_x, reentrant_y, tripolar_n, symmetric
142 
143  if (associated(cs)) call mpp_error(fatal,'Calling oda_init with associated control structure')
144  allocate(cs)
145 ! Use ens1 parameters , this could be changed at a later time
146 ! if it were desirable to have alternate parameters, e.g. for the grid
147 ! for the analysis
148  call get_mom_input(pf,dirs,ensemble_num=0)
149 
150  call unit_scaling_init(pf, cs%US)
151 
152  call get_param(pf, "MOM", "ASSIM_METHOD", assim_method, &
153  "String which determines the data assimilation method "//&
154  "Valid methods are: \'EAKF\',\'OI\', and \'NO_ASSIM\'", default='NO_ASSIM')
155  call get_param(pf, "MOM", "ASSIM_FREQUENCY", cs%assim_frequency, &
156  "data assimilation frequency in hours")
157  call get_param(pf, "MOM", "USE_REGRIDDING", cs%use_ALE_algorithm , &
158  "If True, use the ALE algorithm (regridding/remapping).\n"//&
159  "If False, use the layered isopycnal algorithm.", default=.false. )
160  call get_param(pf, "MOM", "REENTRANT_X", cs%reentrant_x, &
161  "If true, the domain is zonally reentrant.", default=.true.)
162  call get_param(pf, "MOM", "REENTRANT_Y", cs%reentrant_y, &
163  "If true, the domain is meridionally reentrant.", &
164  default=.false.)
165  call get_param(pf,"MOM", "TRIPOLAR_N", cs%tripolar_N, &
166  "Use tripolar connectivity at the northern edge of the "//&
167  "domain. With TRIPOLAR_N, NIGLOBAL must be even.", &
168  default=.false.)
169  call get_param(pf,"MOM", "NIGLOBAL", cs%ni, &
170  "The total number of thickness grid points in the "//&
171  "x-direction in the physical domain.")
172  call get_param(pf,"MOM", "NJGLOBAL", cs%nj, &
173  "The total number of thickness grid points in the "//&
174  "y-direction in the physical domain.")
175  call get_param(pf, 'MOM', "INPUTDIR", inputdir)
176  inputdir = slasher(inputdir)
177 
178  select case(lowercase(trim(assim_method)))
179  case('eakf')
180  cs%assim_method = eakf_assim
181  case('oi')
182  cs%assim_method = oi_assim
183  case('no_assim')
184  cs%assim_method = no_assim
185  case default
186  call mpp_error(fatal,'Invalid assimilation method provided')
187  end select
188 
189  ens_info = get_ensemble_size()
190  cs%ensemble_size = ens_info(1)
191  npes_pm=ens_info(3)
192  cs%ensemble_id = get_ensemble_id()
193  !! Switch to global pelist
194  allocate(cs%ensemble_pelist(cs%ensemble_size,npes_pm))
195  allocate(cs%filter_pelist(cs%ensemble_size*npes_pm))
196  call get_ensemble_pelist(cs%ensemble_pelist,'ocean')
197  call get_ensemble_filter_pelist(cs%filter_pelist,'ocean')
198 
199  call set_current_pelist(cs%filter_pelist)
200 
201  allocate(cs%domains(cs%ensemble_size))
202  cs%domains(cs%ensemble_id)%mpp_domain => g%Domain%mpp_domain
203  do n=1,cs%ensemble_size
204  if (.not. associated(cs%domains(n)%mpp_domain)) allocate(cs%domains(n)%mpp_domain)
205  call set_root_pe(cs%ensemble_pelist(n,1))
206  call mpp_broadcast_domain(cs%domains(n)%mpp_domain)
207  enddo
208  call set_root_pe(cs%filter_pelist(1))
209  allocate(cs%Grid)
210  ! params NIHALO_ODA, NJHALO_ODA set the DA halo size
211  call mom_domains_init(cs%Grid%Domain,pf,param_suffix='_ODA')
212  allocate(hi)
213  call hor_index_init(cs%Grid%Domain, hi, pf, &
214  local_indexing=.false.) ! Use global indexing for DA
215  call verticalgridinit( pf, cs%GV, cs%US )
216  allocate(dg)
217  call create_dyn_horgrid(dg, hi)
218  call clone_mom_domain(cs%Grid%Domain, dg%Domain,symmetric=.false.)
219  call set_grid_metrics(dg,pf)
220  call mom_initialize_topography(dg%bathyT,dg%max_depth,dg,pf)
221  call mom_initialize_coord(cs%GV, cs%US, pf, .false., &
222  dirs%output_directory, tv_dummy, dg%max_depth)
223  call ale_init(pf, cs%GV, cs%US, dg%max_depth, cs%ALE_CS)
224  call mom_grid_init(cs%Grid, pf, global_indexing=.true.)
225  call ale_updateverticalgridtype(cs%ALE_CS, cs%GV)
226  call copy_dyngrid_to_mom_grid(dg, cs%Grid, cs%US)
227  cs%mpp_domain => cs%Grid%Domain%mpp_domain
228  cs%Grid%ke = cs%GV%ke
229  cs%nk = cs%GV%ke
230  ! initialize storage for prior and posterior
231  allocate(cs%Ocean_prior)
232  call init_ocean_ensemble(cs%Ocean_prior,cs%Grid,cs%GV,cs%ensemble_size)
233  allocate(cs%Ocean_posterior)
234  call init_ocean_ensemble(cs%Ocean_posterior,cs%Grid,cs%GV,cs%ensemble_size)
235  allocate(cs%tv)
236 
237  call get_param(pf, 'oda_driver', "REGRIDDING_COORDINATE_MODE", coord_mode, &
238  "Coordinate mode for vertical regridding.", &
239  default="ZSTAR", fail_if_missing=.false.)
240  call initialize_regridding(cs%regridCS, cs%GV, cs%US, dg%max_depth,pf,'oda_driver',coord_mode,'','')
241  call initialize_remapping(cs%remapCS,'PLM')
242  call set_regrid_params(cs%regridCS, min_thickness=0.)
243  call mpp_get_data_domain(g%Domain%mpp_domain,isd,ied,jsd,jed)
244  if (.not. associated(cs%h)) then
245  allocate(cs%h(isd:ied,jsd:jed,cs%GV%ke)); cs%h(:,:,:)=0.0
246  ! assign thicknesses
247  call ale_initthicknesstocoord(cs%ALE_CS,g,cs%GV,cs%h)
248  endif
249  allocate(cs%tv%T(isd:ied,jsd:jed,cs%GV%ke)); cs%tv%T(:,:,:)=0.0
250  allocate(cs%tv%S(isd:ied,jsd:jed,cs%GV%ke)); cs%tv%S(:,:,:)=0.0
251 
252  call set_axes_info(cs%Grid, cs%GV, cs%US, pf, cs%diag_cs, set_vertical=.true.)
253  do n=1,cs%ensemble_size
254  write(fldnam,'(a,i2.2)') 'temp_prior_',n
255  cs%Ocean_prior%id_t(n)=register_diag_field('ODA',trim(fldnam),cs%diag_cs%axesTL%handles,time, &
256  'ocean potential temperature','degC')
257  write(fldnam,'(a,i2.2)') 'salt_prior_',n
258  cs%Ocean_prior%id_s(n)=register_diag_field('ODA',trim(fldnam),cs%diag_cs%axesTL%handles,time, &
259  'ocean salinity','psu')
260  write(fldnam,'(a,i2.2)') 'temp_posterior_',n
261  cs%Ocean_posterior%id_t(n)=register_diag_field('ODA',trim(fldnam),cs%diag_cs%axesTL%handles,time, &
262  'ocean potential temperature','degC')
263  write(fldnam,'(a,i2.2)') 'salt_posterior_',n
264  cs%Ocean_posterior%id_s(n)=register_diag_field('ODA',trim(fldnam),cs%diag_cs%axesTL%handles,time, &
265  'ocean salinity','psu')
266  enddo
267 
268  call mpp_get_data_domain(cs%mpp_domain,isd,ied,jsd,jed)
269  allocate(cs%oda_grid)
270  cs%oda_grid%x => cs%Grid%geolonT
271  cs%oda_grid%y => cs%Grid%geolatT
272 
273  call get_param(pf, 'oda_driver', "BASIN_FILE", basin_file, &
274  "A file in which to find the basin masks, in variable 'basin'.", &
275  default="basin.nc")
276  basin_file = trim(inputdir) // trim(basin_file)
277  allocate(cs%oda_grid%basin_mask(isd:ied,jsd:jed))
278  cs%oda_grid%basin_mask(:,:) = 0.0
279  call mom_read_data(basin_file,'basin',cs%oda_grid%basin_mask,cs%Grid%domain, timelevel=1)
280 
281 ! get global grid information from ocean_model
282  allocate(t_grid)
283  allocate(t_grid%x(cs%ni,cs%nj))
284  allocate(t_grid%y(cs%ni,cs%nj))
285  allocate(t_grid%basin_mask(cs%ni,cs%nj))
286  call mpp_global_field(cs%mpp_domain, cs%Grid%geolonT, t_grid%x)
287  call mpp_global_field(cs%mpp_domain, cs%Grid%geolatT, t_grid%y)
288  call mpp_global_field(cs%mpp_domain, cs%oda_grid%basin_mask, t_grid%basin_mask)
289  t_grid%ni = cs%ni
290  t_grid%nj = cs%nj
291  t_grid%nk = cs%nk
292  allocate(t_grid%mask(cs%ni,cs%nj,cs%nk))
293  allocate(t_grid%z(cs%ni,cs%nj,cs%nk))
294  allocate(global2d(cs%ni,cs%nj))
295  allocate(global2d_old(cs%ni,cs%nj))
296  t_grid%mask(:,:,:) = 0.0
297  t_grid%z(:,:,:) = 0.0
298 
299  do k = 1, cs%nk
300  call mpp_global_field(g%Domain%mpp_domain, cs%h(:,:,k), global2d)
301  do i=1, cs%ni; do j=1, cs%nj
302  if ( global2d(i,j) > 1 ) then
303  t_grid%mask(i,j,k) = 1.0
304  endif
305  enddo ; enddo
306  if (k == 1) then
307  t_grid%z(:,:,k) = global2d/2
308  else
309  t_grid%z(:,:,k) = t_grid%z(:,:,k-1) + (global2d + global2d_old)/2
310  endif
311  global2d_old = global2d
312  enddo
313 
314  call ocean_da_core_init(cs%mpp_domain, t_grid, cs%Profiles, time)
315 
316  cs%Time=time
317  !! switch back to ensemble member pelist
318  call set_current_pelist(cs%ensemble_pelist(cs%ensemble_id,:))
319 end subroutine init_oda
320 
321 !> Copy ensemble member tracers to ensemble vector.
322 subroutine set_prior_tracer(Time, G, GV, h, tv, CS)
323  type(time_type), intent(in) :: time !< The current model time
324  type(ocean_grid_type), pointer :: g !< domain and grid information for ocean model
325  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
326  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
327  type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables
328 
329  type(oda_cs), pointer :: cs !< ocean DA control structure
330  real, dimension(:,:,:), allocatable :: t, s
331  type(ocean_grid_type), pointer :: grid=>null()
332  integer :: i,j, m, n, ss
333  integer :: is, ie, js, je
334  integer :: isc, iec, jsc, jec
335  integer :: isd, ied, jsd, jed
336  integer :: id
337  logical :: used
338 
339  ! return if not time for analysis
340  if (time < cs%Time) return
341 
342  if (.not. associated(cs%Grid)) call mom_error(fatal,'ODA_CS ensemble horizontal grid not associated')
343  if (.not. associated(cs%GV)) call mom_error(fatal,'ODA_CS ensemble vertical grid not associated')
344 
345  !! switch to global pelist
346  call set_current_pelist(cs%filter_pelist)
347  call mom_mesg('Setting prior')
348 
349  isc=cs%Grid%isc;iec=cs%Grid%iec;jsc=cs%Grid%jsc;jec=cs%Grid%jec
350  call mpp_get_compute_domain(cs%domains(cs%ensemble_id)%mpp_domain,is,ie,js,je)
351  call mpp_get_data_domain(cs%domains(cs%ensemble_id)%mpp_domain,isd,ied,jsd,jed)
352  allocate(t(isd:ied,jsd:jed,cs%nk))
353  allocate(s(isd:ied,jsd:jed,cs%nk))
354 
355  do j=js,je; do i=is,ie
356  call remapping_core_h(cs%remapCS, gv%ke, h(i,j,:), tv%T(i,j,:), &
357  cs%nk, cs%h(i,j,:), t(i,j,:))
358  call remapping_core_h(cs%remapCS, gv%ke, h(i,j,:), tv%S(i,j,:), &
359  cs%nk, cs%h(i,j,:), s(i,j,:))
360  enddo ; enddo
361 
362  do m=1,cs%ensemble_size
363  call mpp_redistribute(cs%domains(m)%mpp_domain, t,&
364  cs%mpp_domain, cs%Ocean_prior%T(:,:,:,m), complete=.true.)
365  call mpp_redistribute(cs%domains(m)%mpp_domain, s,&
366  cs%mpp_domain, cs%Ocean_prior%S(:,:,:,m), complete=.true.)
367  if (cs%Ocean_prior%id_t(m)>0) &
368  used=send_data(cs%Ocean_prior%id_t(m), cs%Ocean_prior%T(isc:iec,jsc:jec,:,m), cs%Time)
369  if (cs%Ocean_prior%id_s(m)>0) &
370  used=send_data(cs%Ocean_prior%id_s(m), cs%Ocean_prior%S(isc:iec,jsc:jec,:,m), cs%Time)
371  enddo
372  deallocate(t,s)
373 
374  !! switch back to ensemble member pelist
375  call set_current_pelist(cs%ensemble_pelist(cs%ensemble_id,:))
376 
377  return
378 
379 end subroutine set_prior_tracer
380 
381 !> Returns posterior adjustments or full state
382 !!Note that only those PEs associated with an ensemble member receive data
383 subroutine get_posterior_tracer(Time, CS, h, tv, increment)
384  type(time_type), intent(in) :: time !< the current model time
385  type(oda_cs), pointer :: cs !< ocean DA control structure
386  real, dimension(:,:,:), pointer :: h !< Layer thicknesses [H ~> m or kg m-2]
387  type(thermo_var_ptrs), pointer :: tv !< A structure pointing to various thermodynamic variables
388  logical, optional, intent(in) :: increment !< True if returning increment only
389 
390  type(ocean_control_struct), pointer :: ocean_increment=>null()
391  integer :: i, j, m
392  logical :: used, get_inc
393 
394  ! return if not analysis time (retain pointers for h and tv)
395  if (time < cs%Time) return
396 
397 
398  !! switch to global pelist
399  call set_current_pelist(cs%filter_pelist)
400  call mom_mesg('Getting posterior')
401 
402  get_inc = .true.
403  if (present(increment)) get_inc = increment
404 
405  if (get_inc) then
406  allocate(ocean_increment)
407  call init_ocean_ensemble(ocean_increment,cs%Grid,cs%GV,cs%ensemble_size)
408  ocean_increment%T = cs%Ocean_posterior%T - cs%Ocean_prior%T
409  ocean_increment%S = cs%Ocean_posterior%S - cs%Ocean_prior%S
410  endif
411  do m=1,cs%ensemble_size
412  if (get_inc) then
413  call mpp_redistribute(cs%mpp_domain, ocean_increment%T(:,:,:,m), &
414  cs%domains(m)%mpp_domain, cs%tv%T, complete=.true.)
415  call mpp_redistribute(cs%mpp_domain, ocean_increment%S(:,:,:,m), &
416  cs%domains(m)%mpp_domain, cs%tv%S, complete=.true.)
417  else
418  call mpp_redistribute(cs%mpp_domain, cs%Ocean_posterior%T(:,:,:,m), &
419  cs%domains(m)%mpp_domain, cs%tv%T, complete=.true.)
420  call mpp_redistribute(cs%mpp_domain, cs%Ocean_posterior%S(:,:,:,m), &
421  cs%domains(m)%mpp_domain, cs%tv%S, complete=.true.)
422  endif
423  enddo
424 
425  tv => cs%tv
426  h => cs%h
427  !! switch back to ensemble member pelist
428  call set_current_pelist(cs%ensemble_pelist(cs%ensemble_id,:))
429 
430 end subroutine get_posterior_tracer
431 
432 !> Gather observations and sall ODA routines
433 subroutine oda(Time, CS)
434  type(time_type), intent(in) :: time !< the current model time
435  type(oda_cs), intent(inout) :: cs !< the ocean DA control structure
436 
437  integer :: i, j
438  integer :: m
439  integer :: yr, mon, day, hr, min, sec
440 
441  if ( time >= cs%Time ) then
442 
443  !! switch to global pelist
444  call set_current_pelist(cs%filter_pelist)
445 
446  call get_profiles(time, cs%Profiles, cs%CProfiles)
447 #ifdef ENABLE_ECDA
448  call ensemble_filter(cs%Ocean_prior, cs%Ocean_posterior, cs%CProfiles, cs%kdroot, cs%mpp_domain, cs%oda_grid)
449 #endif
450 
451  !! switch back to ensemble member pelist
452  call set_current_pelist(cs%ensemble_pelist(cs%ensemble_id,:))
453 
454  endif
455 
456  return
457 end subroutine oda
458 
459 !> Finalize DA module
460 subroutine oda_end(CS)
461  type(oda_cs), intent(inout) :: cs !< the ocean DA control structure
462 
463 end subroutine oda_end
464 
465 !> Initialize DA module
466 subroutine init_ocean_ensemble(CS,Grid,GV,ens_size)
467  type(ocean_control_struct), pointer :: CS !< Pointer to ODA control structure
468  type(ocean_grid_type), pointer :: Grid !< Pointer to ocean analysis grid
469  type(verticalgrid_type), pointer :: GV !< Pointer to DA vertical grid
470  integer, intent(in) :: ens_size !< ensemble size
471 
472  integer :: n,is,ie,js,je,nk
473 
474  nk=gv%ke
475  is=grid%isd;ie=grid%ied
476  js=grid%jsd;je=grid%jed
477  cs%ensemble_size=ens_size
478  allocate(cs%T(is:ie,js:je,nk,ens_size))
479  allocate(cs%S(is:ie,js:je,nk,ens_size))
480  allocate(cs%SSH(is:ie,js:je,ens_size))
481  allocate(cs%id_t(ens_size));cs%id_t(:)=-1
482  allocate(cs%id_s(ens_size));cs%id_s(:)=-1
483 ! allocate(CS%U(is:ie,js:je,nk,ens_size))
484 ! allocate(CS%V(is:ie,js:je,nk,ens_size))
485 ! allocate(CS%id_u(ens_size));CS%id_u(:)=-1
486 ! allocate(CS%id_v(ens_size));CS%id_v(:)=-1
487  allocate(cs%id_ssh(ens_size));cs%id_ssh(:)=-1
488 
489  return
490 end subroutine init_ocean_ensemble
491 
492 !> Set the next analysis time
493 subroutine set_analysis_time(Time,CS)
494  type(time_type), intent(in) :: time !< the current model time
495  type(oda_cs), pointer, intent(inout) :: cs !< the DA control structure
496 
497  character(len=160) :: mesg ! The text of an error message
498  integer :: yr, mon, day, hr, min, sec
499 
500  if (time >= cs%Time) then
501  cs%Time=increment_time(cs%Time,cs%assim_frequency*3600)
502 
503  call get_date(time, yr, mon, day, hr, min, sec)
504  write(mesg,*) 'Model Time: ', yr, mon, day, hr, min, sec
505  call mom_mesg("set_analysis_time: "//trim(mesg))
506  call get_date(cs%time, yr, mon, day, hr, min, sec)
507  write(mesg,*) 'Assimilation Time: ', yr, mon, day, hr, min, sec
508  call mom_mesg("set_analysis_time: "//trim(mesg))
509  endif
510  if (cs%Time < time) then
511  call mom_error(fatal, " set_analysis_time: " // &
512  "assimilation interval appears to be shorter than " // &
513  "the model timestep")
514  endif
515  return
516 
517 end subroutine set_analysis_time
518 
519 !> Write observation differences to a file
520 subroutine save_obs_diff(filename,CS)
521  character(len=*), intent(in) :: filename !< name of output file
522  type(oda_cs), pointer, intent(in) :: cs !< pointer to DA control structure
523 
524  integer :: fid ! profile file handle
525  type(ocean_profile_type), pointer :: prof=>null()
526 
527  fid = open_profile_file(trim(filename), nvar=2, thread=mpp_single, fset=mpp_single)
528  prof=>cs%CProfiles
529 
530  !! switch to global pelist
531  !call set_current_pelist(CS%filter_pelist)
532 
533  do while (associated(prof))
534  call write_profile(fid,prof)
535  prof=>prof%cnext
536  enddo
537  call close_profile_file(fid)
538 
539  !! switch back to ensemble member pelist
540  !call set_current_pelist(CS%ensemble_pelist(CS%ensemble_id,:))
541 
542  return
543 end subroutine save_obs_diff
544 
545 
546 !> Apply increments to tracers
547 subroutine apply_oda_tracer_increments(dt,G,tv,h,CS)
548  real, intent(in) :: dt !< The tracer timestep [s]
549  type(ocean_grid_type), intent(in) :: g !< ocean grid structure
550  type(thermo_var_ptrs), intent(inout) :: tv !< A structure pointing to various thermodynamic variables
551  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
552  intent(in) :: h !< layer thickness [H ~> m or kg m-2]
553  type(oda_cs), intent(inout) :: cs !< the data assimilation structure
554 
555 end subroutine apply_oda_tracer_increments
556 
557 !> \namespace MOM_oda_driver_mod
558 !!
559 !! \section section_ODA The Ocean data assimilation (DA) and Ensemble Framework
560 !!
561 !! The DA framework implements ensemble capability in MOM6. Currently, this framework
562 !! is enabled using the cpp directive ENSEMBLE_OCEAN. The ensembles need to be generated
563 !! at the level of the calling routine for oda_init or above. The ensemble instances may
564 !! exist on overlapping or non-overlapping processors. The ensemble information is accessed
565 !! via the FMS ensemble manager. An independent PE layout is used to gather (prior) ensemble
566 !! member information where this information is stored in the ODA control structure. This
567 !! module was developed in collaboration with Feiyu Lu and Tony Rosati in the GFDL prediction
568 !! group for use in their coupled ensemble framework. These interfaces should be suitable for
569 !! interfacing MOM6 to other data assimilation packages as well.
570 
571 end module mom_oda_driver_mod
mom_fixed_initialization::mom_initialize_fixed
subroutine, public mom_initialize_fixed(G, US, OBC, PF, write_geom, output_dir)
MOM_initialize_fixed sets up time-invariant quantities related to MOM6's horizontal grid,...
Definition: MOM_fixed_initialization.F90:56
mom_grid_initialize::set_grid_metrics
subroutine, public set_grid_metrics(G, param_file, US)
set_grid_metrics is used to set the primary values in the model's horizontal grid....
Definition: MOM_grid_initialize.F90:63
mom_regridding::set_regrid_params
subroutine, public set_regrid_params(CS, boundary_extrapolation, min_thickness, old_grid_weight, interp_scheme, depth_of_time_filter_shallow, depth_of_time_filter_deep, compress_fraction, dz_min_surface, nz_fixed_surface, Rho_ML_avg_depth, nlay_ML_to_interior, fix_haloclines, halocline_filt_len, halocline_strat_tol, integrate_downward_for_e, adaptTimeRatio, adaptZoom, adaptZoomCoeff, adaptBuoyCoeff, adaptAlpha, adaptDoMin)
Can be used to set any of the parameters for MOM_regridding.
Definition: MOM_regridding.F90:2218
mom_regridding::initialize_regridding
subroutine, public initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_mode, param_prefix, param_suffix)
Initialization and configures a regridding control structure based on customizable run-time parameter...
Definition: MOM_regridding.F90:177
mom_regridding::regridding_cs
Regridding control structure.
Definition: MOM_regridding.F90:45
mom_remapping::remapping_core_h
subroutine, public remapping_core_h(CS, n0, h0, u0, n1, h1, u1, h_neglect, h_neglect_edge)
Remaps column of values u0 on grid h0 to grid h1 assuming the top edge is aligned.
Definition: MOM_remapping.F90:188
mom_verticalgrid
Provides a transparent vertical ocean grid type and supporting routines.
Definition: MOM_verticalGrid.F90:2
mom_dyn_horgrid::destroy_dyn_horgrid
subroutine, public destroy_dyn_horgrid(G)
Release memory used by the dyn_horgrid_type and related structures.
Definition: MOM_dyn_horgrid.F90:377
mom_oda_driver_mod::oda
subroutine, public oda(Time, CS)
Gather observations and sall ODA routines.
Definition: MOM_oda_driver.F90:434
mom_oda_driver_mod::set_prior_tracer
subroutine, public set_prior_tracer(Time, G, GV, h, tv, CS)
Copy ensemble member tracers to ensemble vector.
Definition: MOM_oda_driver.F90:323
mom_oda_driver_mod::eakf_assim
integer, parameter eakf_assim
DA parameters.
Definition: MOM_oda_driver.F90:108
mom_oda_driver_mod
Interfaces for MOM6 ensembles and data assimilation.
Definition: MOM_oda_driver.F90:2
mom_diag_mediator
The subroutines here provide convenient wrappers to the fms diag_manager interfaces with additional d...
Definition: MOM_diag_mediator.F90:3
mom_variables::thermo_var_ptrs
Pointers to an assortment of thermodynamic fields that may be available, including potential temperat...
Definition: MOM_variables.F90:78
mom_get_input::directories
Container for paths and parameter file names.
Definition: MOM_get_input.F90:20
mom_dyn_horgrid
Contains a shareable dynamic type for describing horizontal grids and metric data and utilty routines...
Definition: MOM_dyn_horgrid.F90:3
mom_oda_driver_mod::init_ocean_ensemble
subroutine init_ocean_ensemble(CS, Grid, GV, ens_size)
Initialize DA module.
Definition: MOM_oda_driver.F90:467
mom_string_functions
Handy functions for manipulating strings.
Definition: MOM_string_functions.F90:2
mom_hor_index::hor_index_init
subroutine, public hor_index_init(Domain, HI, param_file, local_indexing, index_offset)
Sets various index values in a hor_index_type.
Definition: MOM_hor_index.F90:58
mom_coord_initialization::mom_initialize_coord
subroutine, public mom_initialize_coord(GV, US, PF, write_geom, output_dir, tv, max_depth)
MOM_initialize_coord sets up time-invariant quantities related to MOM6's vertical coordinate.
Definition: MOM_coord_initialization.F90:40
mom_diag_mediator::set_axes_info
subroutine, public set_axes_info(G, GV, US, param_file, diag_cs, set_vertical)
Sets up diagnostics axes.
Definition: MOM_diag_mediator.F90:344
mom_regridding::regridding_main
subroutine, public regridding_main(remapCS, CS, G, GV, h, tv, h_new, dzInterface, frac_shelf_h, conv_adjust)
Dispatching regridding routine for orchestrating regridding & remapping.
Definition: MOM_regridding.F90:791
mom_error_handler::mom_mesg
subroutine, public mom_mesg(message, verb, all_print)
This provides a convenient interface for writing an informative comment.
Definition: MOM_error_handler.F90:53
mom_oda_driver_mod::ptr_mpp_domain
A structure with a pointer to a domain2d, to allow for the creation of arrays of pointers.
Definition: MOM_oda_driver.F90:103
mom_remapping::remapping_cs
Container for remapping parameters.
Definition: MOM_remapping.F90:22
mom_file_parser::param_file_type
A structure that can be parsed to read and document run-time parameters.
Definition: MOM_file_parser.F90:54
mom_file_parser::get_param
An overloaded interface to read and log the values of various types of parameters.
Definition: MOM_file_parser.F90:102
mom_hor_index
Defines the horizontal index type (hor_index_type) used for providing index ranges.
Definition: MOM_hor_index.F90:2
mom_oda_driver_mod::no_assim
integer, parameter no_assim
DA parameters.
Definition: MOM_oda_driver.F90:108
mom_oda_driver_mod::oda_end
subroutine, public oda_end(CS)
Finalize DA module.
Definition: MOM_oda_driver.F90:461
mom_io
This module contains I/O framework code.
Definition: MOM_io.F90:2
mom_fixed_initialization::mom_initialize_topography
subroutine, public mom_initialize_topography(D, max_depth, G, PF, US)
MOM_initialize_topography makes the appropriate call to set up the bathymetry. At this point the topo...
Definition: MOM_fixed_initialization.F90:174
mom_transcribe_grid::copy_mom_grid_to_dyngrid
subroutine, public copy_mom_grid_to_dyngrid(oG, dG, US)
Copies information from an ocean_grid_type into a dynamic (shared) horizontal grid type.
Definition: MOM_transcribe_grid.F90:167
mom_get_input
Reads the only Fortran name list needed to boot-strap the model.
Definition: MOM_get_input.F90:6
mom_oda_driver_mod::set_analysis_time
subroutine, public set_analysis_time(Time, CS)
Set the next analysis time.
Definition: MOM_oda_driver.F90:494
mom_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
mom_verticalgrid::verticalgridinit
subroutine, public verticalgridinit(param_file, GV, US)
Allocates and initializes the ocean model vertical grid structure.
Definition: MOM_verticalGrid.F90:76
mom_remapping
Provides column-wise vertical remapping functions.
Definition: MOM_remapping.F90:2
mom_transcribe_grid
Module with routines for copying information from a shared dynamic horizontal grid to an ocean-specif...
Definition: MOM_transcribe_grid.F90:3
mom_domains::clone_mom_domain
Copy one MOM_domain_type into another.
Definition: MOM_domains.F90:94
mom_ale::ale_initthicknesstocoord
subroutine, public ale_initthicknesstocoord(CS, G, GV, h)
Set h to coordinate values for fixed coordinate systems.
Definition: MOM_ALE.F90:1286
mom_grid::mom_grid_init
subroutine, public mom_grid_init(G, param_file, US, HI, global_indexing, bathymetry_at_vel)
MOM_grid_init initializes the ocean grid array sizes and grid memory.
Definition: MOM_grid.F90:184
mom_ale::ale_updateverticalgridtype
subroutine, public ale_updateverticalgridtype(CS, GV)
Update the vertical grid type with ALE information. This subroutine sets information in the verticalG...
Definition: MOM_ALE.F90:1235
mom_coord_initialization
Initializes fixed aspects of the related to its vertical coordinate.
Definition: MOM_coord_initialization.F90:2
mom_verticalgrid::verticalgrid_type
Describes the vertical ocean grid, including unit conversion factors.
Definition: MOM_verticalGrid.F90:24
mom_domains
Describes the decomposed MOM domain and has routines for communications across PEs.
Definition: MOM_domains.F90:2
mom_dyn_horgrid::create_dyn_horgrid
subroutine, public create_dyn_horgrid(G, HI, bathymetry_at_vel)
Allocate memory used by the dyn_horgrid_type and related structures.
Definition: MOM_dyn_horgrid.F90:176
mom_ale
This module contains the main regridding routines.
Definition: MOM_ALE.F90:9
mom_variables
Provides transparent structures with groups of MOM6 variables and supporting routines.
Definition: MOM_variables.F90:2
mom_ale::ale_init
subroutine, public ale_init(param_file, GV, US, max_depth, CS)
This routine is typically called (from initialize_MOM in file MOM.F90) before the main time integrati...
Definition: MOM_ALE.F90:141
mom_oda_driver_mod::save_obs_diff
subroutine, public save_obs_diff(filename, CS)
Write observation differences to a file.
Definition: MOM_oda_driver.F90:521
mom_oda_driver_mod::init_oda
subroutine, public init_oda(Time, G, GV, CS)
initialize First_guess (prior) and Analysis grid information for all ensemble members
Definition: MOM_oda_driver.F90:116
mom_io::mom_read_data
Read a data field from a file.
Definition: MOM_io.F90:74
mom_fixed_initialization
Initializes fixed aspects of the model, such as horizontal grid metrics, topography and Coriolis.
Definition: MOM_fixed_initialization.F90:3
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
mom_ale::ale_cs
ALE control structure.
Definition: MOM_ALE.F90:65
mom_oda_driver_mod::apply_oda_tracer_increments
subroutine, public apply_oda_tracer_increments(dt, G, tv, h, CS)
Apply increments to tracers.
Definition: MOM_oda_driver.F90:548
mom_get_input::get_mom_input
subroutine, public get_mom_input(param_file, dirs, check_params, default_input_filename, ensemble_num)
Get the names of the I/O directories and initialization file. Also calls the subroutine that opens ru...
Definition: MOM_get_input.F90:34
mom_hor_index::hor_index_type
Container for horizontal index ranges for data, computational and global domains.
Definition: MOM_hor_index.F90:15
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.F90:2
mom_unit_scaling
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Definition: MOM_unit_scaling.F90:2
mom_error_handler::is_root_pe
logical function, public is_root_pe()
This returns .true. if the current PE is the root PE.
Definition: MOM_error_handler.F90:44
mom_domains::mom_domains_init
subroutine, public mom_domains_init(MOM_dom, param_file, symmetric, static_memory, NIHALO, NJHALO, NIGLOBAL, NJGLOBAL, NIPROC, NJPROC, min_halo, domain_name, include_name, param_suffix)
MOM_domains_init initalizes a MOM_domain_type variable, based on the information read in from a param...
Definition: MOM_domains.F90:1155
mom_grid_initialize
Initializes horizontal grid.
Definition: MOM_grid_initialize.F90:2
mom_transcribe_grid::copy_dyngrid_to_mom_grid
subroutine, public copy_dyngrid_to_mom_grid(dG, oG, US)
Copies information from a dynamic (shared) horizontal grid type into an ocean_grid_type.
Definition: MOM_transcribe_grid.F90:23
mom_error_handler::mom_error
subroutine, public mom_error(level, message, all_print)
This provides a convenient interface for writing an mpp_error message with run-time filter based on a...
Definition: MOM_error_handler.F90:72
mom_oda_driver_mod::oi_assim
integer, parameter oi_assim
DA parameters.
Definition: MOM_oda_driver.F90:108
mom_oda_driver_mod::get_posterior_tracer
subroutine, public get_posterior_tracer(Time, CS, h, tv, increment)
Returns posterior adjustments or full state Note that only those PEs associated with an ensemble memb...
Definition: MOM_oda_driver.F90:384
mom_string_functions::lowercase
character(len=len(input_string)) function, public lowercase(input_string)
Return a string in which all uppercase letters have been replaced by their lowercase counterparts.
Definition: MOM_string_functions.F90:24
mom_regridding
Generates vertical grids as part of the ALE algorithm.
Definition: MOM_regridding.F90:2
mom_domains::mom_domain_type
The MOM_domain_type contains information about the domain decompositoin.
Definition: MOM_domains.F90:99
mom_remapping::initialize_remapping
subroutine, public initialize_remapping(CS, remapping_scheme, boundary_extrapolation, check_reconstruction, check_remapping, force_bounds_in_subcell, answers_2018)
Constructor for remapping control structure.
Definition: MOM_remapping.F90:1547
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
mom_dyn_horgrid::dyn_horgrid_type
Describes the horizontal ocean grid with only dynamic memory arrays.
Definition: MOM_dyn_horgrid.F90:23
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:26
mom_unit_scaling::unit_scaling_init
subroutine, public unit_scaling_init(param_file, US)
Allocates and initializes the ocean model unit scaling type.
Definition: MOM_unit_scaling.F90:44
mom_diag_mediator::diag_ctrl
The following data type a list of diagnostic fields an their variants, as well as variables that cont...
Definition: MOM_diag_mediator.F90:239
mom_oda_driver_mod::oda_cs
Control structure that contains a transpose of the ocean state across ensemble members.
Definition: MOM_oda_driver.F90:62
mom_file_parser::read_param
An overloaded interface to read various types of parameters.
Definition: MOM_file_parser.F90:90