initialize First_guess (prior) and Analysis grid information for all ensemble members
117 type(time_type),
intent(in) :: Time
118 type(ocean_grid_type),
pointer :: G
119 type(verticalGrid_type),
intent(in) :: GV
120 type(ODA_CS),
pointer,
intent(inout) :: CS
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
128 type(grid_type),
pointer :: T_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
143 if (
associated(cs))
call mpp_error(fatal,
'Calling oda_init with associated control structure')
148 call get_mom_input(pf,dirs,ensemble_num=0)
150 call unit_scaling_init(pf, cs%US)
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.", &
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.", &
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)
178 select case(lowercase(trim(assim_method)))
180 cs%assim_method = eakf_assim
182 cs%assim_method = oi_assim
184 cs%assim_method = no_assim
186 call mpp_error(fatal,
'Invalid assimilation method provided')
189 ens_info = get_ensemble_size()
190 cs%ensemble_size = ens_info(1)
192 cs%ensemble_id = get_ensemble_id()
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')
199 call set_current_pelist(cs%filter_pelist)
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)
208 call set_root_pe(cs%filter_pelist(1))
211 call mom_domains_init(cs%Grid%Domain,pf,param_suffix=
'_ODA')
213 call hor_index_init(cs%Grid%Domain, hi, pf, &
214 local_indexing=.false.)
215 call verticalgridinit( pf, cs%GV, cs%US )
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)
227 cs%mpp_domain => cs%Grid%Domain%mpp_domain
228 cs%Grid%ke = cs%GV%ke
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)
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
247 call ale_initthicknesstocoord(cs%ALE_CS,g,cs%GV,cs%h)
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
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')
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
273 call get_param(pf,
'oda_driver',
"BASIN_FILE", basin_file, &
274 "A file in which to find the basin masks, in variable 'basin'.", &
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)
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)
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
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
307 t_grid%z(:,:,k) = global2d/2
309 t_grid%z(:,:,k) = t_grid%z(:,:,k-1) + (global2d + global2d_old)/2
311 global2d_old = global2d
314 call ocean_da_core_init(cs%mpp_domain, t_grid, cs%Profiles, time)
318 call set_current_pelist(cs%ensemble_pelist(cs%ensemble_id,:))