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
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
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
54 implicit none ;
private
59 #include <MOM_memory.h>
63 type(ocean_control_struct),
pointer :: ocean_prior=> null()
64 type(ocean_control_struct),
pointer :: ocean_posterior=> null()
74 type(domain2d),
pointer :: mpp_domain => null()
75 type(grid_type),
pointer :: oda_grid
76 real,
pointer,
dimension(:,:,:) :: h => null()
80 logical :: reentrant_x
81 logical :: reentrant_y
84 integer :: assim_method
85 integer :: ensemble_size
86 integer :: ensemble_id = 0
87 integer,
pointer,
dimension(:,:) :: ensemble_pelist
88 integer,
pointer,
dimension(:) :: filter_pelist
89 integer :: assim_frequency
91 type(ocean_profile_type),
pointer :: profiles => null()
92 type(ocean_profile_type),
pointer :: cprofiles => null()
93 type(kd_root),
pointer :: kdroot => null()
95 logical :: use_ale_algorithm
98 type(time_type) :: time
104 type(domain2d),
pointer :: mpp_domain => null()
115 subroutine init_oda(Time, G, GV, CS)
117 type(time_type),
intent(in) :: time
120 type(
oda_cs),
pointer,
intent(inout) :: cs
128 type(grid_type),
pointer :: t_grid
129 real,
dimension(:,:),
allocatable :: global2d, global2d_old
130 real,
dimension(:),
allocatable :: lon1d, lat1d, glon1d, glat1d
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')
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)))
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))
214 local_indexing=.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)
227 cs%mpp_domain => cs%Grid%Domain%mpp_domain
228 cs%Grid%ke = cs%GV%ke
231 allocate(cs%Ocean_prior)
233 allocate(cs%Ocean_posterior)
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,
'',
'')
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
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,:))
323 type(time_type),
intent(in) :: time
326 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
329 type(
oda_cs),
pointer :: cs
330 real,
dimension(:,:,:),
allocatable :: t, s
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
340 if (time < cs%Time)
return
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')
346 call set_current_pelist(cs%filter_pelist)
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))
355 do j=js,je;
do i=is,ie
357 cs%nk, cs%h(i,j,:), t(i,j,:))
359 cs%nk, cs%h(i,j,:), s(i,j,:))
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)
375 call set_current_pelist(cs%ensemble_pelist(cs%ensemble_id,:))
384 type(time_type),
intent(in) :: time
385 type(
oda_cs),
pointer :: cs
386 real,
dimension(:,:,:),
pointer :: h
388 logical,
optional,
intent(in) :: increment
390 type(ocean_control_struct),
pointer :: ocean_increment=>null()
392 logical :: used, get_inc
395 if (time < cs%Time)
return
399 call set_current_pelist(cs%filter_pelist)
403 if (
present(increment)) get_inc = increment
406 allocate(ocean_increment)
408 ocean_increment%T = cs%Ocean_posterior%T - cs%Ocean_prior%T
409 ocean_increment%S = cs%Ocean_posterior%S - cs%Ocean_prior%S
411 do m=1,cs%ensemble_size
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.)
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.)
428 call set_current_pelist(cs%ensemble_pelist(cs%ensemble_id,:))
433 subroutine oda(Time, CS)
434 type(time_type),
intent(in) :: time
435 type(
oda_cs),
intent(inout) :: cs
439 integer :: yr, mon, day, hr, min, sec
441 if ( time >= cs%Time )
then
444 call set_current_pelist(cs%filter_pelist)
446 call get_profiles(time, cs%Profiles, cs%CProfiles)
448 call ensemble_filter(cs%Ocean_prior, cs%Ocean_posterior, cs%CProfiles, cs%kdroot, cs%mpp_domain, cs%oda_grid)
452 call set_current_pelist(cs%ensemble_pelist(cs%ensemble_id,:))
467 type(ocean_control_struct),
pointer :: CS
470 integer,
intent(in) :: ens_size
472 integer :: n,is,ie,js,je,nk
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
487 allocate(cs%id_ssh(ens_size));cs%id_ssh(:)=-1
494 type(time_type),
intent(in) :: time
495 type(
oda_cs),
pointer,
intent(inout) :: cs
497 character(len=160) :: mesg
498 integer :: yr, mon, day, hr, min, sec
500 if (time >= cs%Time)
then
501 cs%Time=increment_time(cs%Time,cs%assim_frequency*3600)
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))
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")
521 character(len=*),
intent(in) :: filename
522 type(
oda_cs),
pointer,
intent(in) :: cs
525 type(ocean_profile_type),
pointer :: prof=>null()
527 fid = open_profile_file(trim(filename), nvar=2, thread=mpp_single, fset=mpp_single)
533 do while (
associated(prof))
534 call write_profile(fid,prof)
537 call close_profile_file(fid)
548 real,
intent(in) :: dt
551 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
553 type(
oda_cs),
intent(inout) :: cs