da_solve.f90
References to this file elsewhere.
1 subroutine da_solve ( grid , config_flags)
2
3 !-----------------------------------------------------------------------
4 ! Purpose: TBD
5 !-----------------------------------------------------------------------
6
7 use module_configure, only : grid_config_rec_type
8 use module_domain, only : domain
9 use module_driver_constants, only : max_comms
10 use module_symbols_util, only : wrfu_finalize
11
12 use module_radiance, only : satinfo, time_slots
13 #ifdef RTTOV
14 use module_radiance, only : coefs, sensor_descriptor
15 #endif
16
17 use module_state_description, only : num_moist, num_a_moist, num_g_moist, &
18 num_scalar, num_a_scalar, num_g_scalar
19 use module_tiles, only : set_tiles
20
21 #ifdef DM_PARALLEL
22 use module_dm, only : invalid_message_value, glen,setup_xpose_rsl, &
23 add_msg_24pt_real,reset_msgs_24pt,stencil_24pt,setup_halo_rsl, &
24 reset_msgs_xpose, add_msg_xpose_real, define_xpose, add_msg_24pt_integer
25 #endif
26
27 use da_control, only : trace_use, comm, ierr, ids,ide,jds,jde,kds,kde, &
28 ips,ipe, jps,jpe, vert_corr, sin_xle, testing_wrfvar, use_rad, &
29 calc_w_increment, var4d_coupling_disk_simul, var4d_coupling, &
30 write_oa_rad_ascii, var4d, cos_xls, vertical_ip, use_radarobs, stdout, &
31 sin_xls, rf_passes, ntmax, rootproc,test_transforms,global, &
32 cos_xle,anal_type_qcobs,check_max_iv,anal_type_randomcv,cv_options_hum, &
33 max_ext_its,anal_type_verify, start_x, start_y,coarse_ix, coarse_jy, &
34 rtm_option, rtm_option_crtm, rtm_option_rttov,read_biascoef
35 use da_define_structures, only : y_type, j_type, ob_type, be_type, &
36 xbx_type,da_deallocate_background_errors,da_initialize_cv, &
37 da_zero_vp_type,da_allocate_y,da_deallocate_observations, &
38 da_deallocate_y
39 use da_minimisation, only : da_get_innov_vector,da_minimise_cg, &
40 da_write_diagnostics
41 use da_radiance1, only : da_write_oa_rad_ascii
42 use da_obs_io, only : da_write_filtered_obs
43 use da_par_util, only : da_system,da_copy_tile_dims,da_copy_dims
44 use da_physics, only : da_uvprho_to_w_lin
45 use da_reporting, only : message, da_warning
46 use da_setup_structures, only : da_setup_obs_structures, &
47 da_setup_background_errors,da_setup_flow_predictors
48 use da_test, only : da_check
49 use da_tools1, only : da_get_unit, da_free_unit
50 use da_tracing, only : da_trace_entry, da_trace_exit, da_trace
51 use da_transfer_model, only : da_transfer_xatoanalysis,da_setup_firstguess
52 use da_vtox_transforms, only : da_transform_vtox
53 use da_wrf_interfaces, only : wrf_error_fatal3
54
55 implicit none
56
57 type (domain), intent(inout) :: grid
58 type (grid_config_rec_type), intent(inout) :: config_flags
59
60 type (xbx_type) :: xbx ! For header & non-grid arrays.
61 type (be_type) :: be ! Background error structure.
62 real, allocatable :: cvt(:) ! Control variable structure.
63 real, allocatable :: xhat(:) ! Control variable structure.
64 type (y_type) :: ob ! Observation structure.
65 type (ob_type) :: iv ! Obs. increment structure.
66 type (y_type) :: re ! Residual (o-a) structure.
67 type (y_type) :: y ! y = H(x_inc) structure.
68 integer :: it ! External loop counter.
69 type (j_type) :: j ! Cost function.
70
71 integer :: cv_size, i
72 real :: j_grad_norm_target ! Target j norm.
73
74 #ifdef DM_PARALLEL
75 integer :: wrf_done_unit
76 #endif
77
78 if (trace_use) call da_trace_entry("da_solve")
79
80 #ifdef DM_PARALLEL
81 call mpi_barrier(comm,ierr)
82 #endif
83
84 !---------------------------------------------------------------------------
85 ! If it is verification run set check_max_iv as .false.
86 !---------------------------------------------------------------------------
87
88 if (anal_type_verify) then
89 check_max_iv = .false.
90 ntmax=0
91 end if
92
93 if (cv_options_hum < 1 .or. cv_options_hum > 3) then
94 write(unit=message(1),fmt='(A,I3)') &
95 'Invalid cv_options_hum = ', cv_options_hum
96 call wrf_error_fatal(message(1:1))
97 end if
98
99 if (vert_corr == 2) then
100 if (vertical_ip < 0 .or. vertical_ip > 2) then
101 write (unit=message(1),fmt='(A,I3)') &
102 'Invalid vertical_ip = ', &
103 vertical_ip
104 call da_warning(__FILE__,__LINE__,message(1:1))
105 end if
106 end if
107
108 if (0.5 * real(rf_passes) /= real(rf_passes / 2)) then
109 write(unit=stdout,fmt='(A,I4,A)')'rf_passes = ', &
110 rf_passes, ' .Should be even.'
111 rf_passes = int(real(rf_passes / 2))
112 write(unit=stdout,fmt='(A,I4)') 'Resetting rf_passes = ', rf_passes
113 end if
114
115 if (anal_type_randomcv) then
116 ntmax = 0
117 write(unit=stdout,fmt='(a)')' Resetting ntmax = 0 for analysis_type = randomcv'
118 end if
119
120 !---------------------------------------------------------------------------
121 ! [2.0] Initialise wrfvar parameters:
122 !---------------------------------------------------------------------------
123
124 call da_solve_init(grid &
125 #include "em_actual_new_args.inc"
126 )
127
128 !---------------------------------------------------------------------------
129 ! [3.0] Set up first guess field (grid%xb):
130 !---------------------------------------------------------------------------
131
132 call da_setup_firstguess( xbx, grid)
133
134 !---------------------------------------------------------------------------
135 ! [4.0] Set up observations (ob):
136 !---------------------------------------------------------------------------
137
138 call da_setup_obs_structures( grid, ob, iv)
139
140 if (use_rad) then
141 allocate (j % jo % rad(1:iv%num_inst))
142 do i=1,iv%num_inst
143 allocate (j % jo % rad(i) % jo_ichan(iv%instid(i)%nchan))
144 allocate (j % jo % rad(i) % num_ichan(iv%instid(i)%nchan))
145 end do
146 end if
147
148 !---------------------------------------------------------------------------
149 ! [5.0] Set up background errors (be):
150 !---------------------------------------------------------------------------
151
152 call da_setup_background_errors( grid%xb, be)
153 cv_size = be % cv % size
154
155 !---------------------------------------------------------------------------
156 ! [6.0] Set up ensemble perturbation input:
157 !---------------------------------------------------------------------------
158
159 grid%ep % ne = be % ne
160 if (be % ne > 0) then
161 call da_setup_flow_predictors( ide, jde, kde, be % ne, grid%ep )
162 end if
163
164 !---------------------------------------------------------------------------
165 ! [7.0] Setup control variable (cv):
166 !---------------------------------------------------------------------------
167
168 allocate (cvt(1:cv_size))
169 allocate (xhat(1:cv_size))
170 call da_initialize_cv (cv_size, cvt)
171 call da_initialize_cv (cv_size, xhat)
172
173 call da_zero_vp_type (grid%vv)
174 call da_zero_vp_type (grid%vp)
175
176 !---------------------------------------------------------------------------
177 ! [8] Outerloop
178 !---------------------------------------------------------------------------
179
180 j_grad_norm_target = 1.0
181 do it = 1, max_ext_its
182
183 ! [8.1] Calculate nonlinear model trajectory
184
185 if (var4d) then
186 call da_trace ("da_solve","Starting da_run_wrf_nl.ksh")
187 #ifdef DM_PARALLEL
188 if (var4d_coupling == var4d_coupling_disk_simul) then
189 call da_system ("da_run_wrf_nl.ksh pre")
190 ! call da_system("./wrf.exe -rmpool 1")
191 if (rootproc) then
192 call da_system("rm -rf wrfnl_done")
193 call da_system("touch wrfnl_go_ahead")
194 call da_get_unit(wrf_done_unit)
195 do while ( .true. )
196 open(wrf_done_unit,file="wrfnl_done",status="old",err=303)
197 close(wrf_done_unit)
198 exit
199 303 continue
200 call da_system("sync")
201 call da_system("sleep 1")
202 end do
203 call da_free_unit(wrf_done_unit)
204 end if
205 ! Wait until PE thinks NL model has finished
206 call mpi_barrier( comm, ierr )
207 call da_system("da_run_wrf_nl.ksh post")
208 call da_system("touch wrfnl_stop_now")
209 end if
210 #else
211 call da_system("da_run_wrf_nl.ksh")
212 #endif
213 call da_trace("da_solve","Finished da_run_wrf_nl.ksh")
214 end if
215
216 ! [8.2] Calculate innovation vector (O-B):
217
218 call da_get_innov_vector( it, ob, iv, grid , config_flags)
219
220 if (test_transforms) then
221 call da_check( cv_size, grid%xb, xbx, be, grid%ep, iv, &
222 grid%xa, grid%vv, grid%vp, grid%xp, y)
223 call wrfu_finalize
224 call wrf_shutdown
225 stop
226 end if
227
228 if (testing_wrfvar) then
229 call da_check( cv_size, grid%xb, xbx, be, grid%ep, iv, &
230 grid%xa, grid%vv, grid%vp, grid%xp, y)
231 call wrfu_finalize
232 call wrf_shutdown
233 stop
234 end if
235
236 ! [8.3] Interpolate x_g to low resolution grid
237
238 ! [8.4] Minimize cost function:
239
240 call da_allocate_y( iv, re )
241 call da_allocate_y( iv, y )
242
243 call da_minimise_cg( grid, config_flags, &
244 it, be % cv % size, &
245 xbx, be, iv, &
246 j_grad_norm_target, xhat, cvt, &
247 re, y, j)
248
249 !------------------------------------------------------------------------
250
251 ! [8.5] Update latest analysis solution:
252
253 call da_transform_vtox( cv_size, grid%xb, xbx, be, grid%ep, xhat, &
254 grid%vv, grid%vp, grid%xp, grid%xa)
255
256 ! [8.6] Only when use_RadarObs = .false. and calc_w_increment =.true.,
257 ! the w_increment need to be diagnosed:
258
259 if (calc_w_increment .and. .not. use_RadarObs) then
260 call da_uvprho_to_w_lin( grid%xb, grid%xa, grid%xp)
261
262 call wrf_dm_halo(grid%xp%domdesc,grid%xp%comms,grid%xp%halo_id13)
263 end if
264
265 ! [8.7] Write out diagnostics
266
267 call da_write_diagnostics( ob, iv, re, y, grid%xp, grid%xa, j )
268
269
270 ! Write "clean" QCed observations if requested:
271 if (anal_type_qcobs) then
272 if (it == 1) then
273 call da_write_filtered_obs(ob, iv, grid%xb, grid%xp, &
274 grid%moad_cen_lat, grid%stand_lon,&
275 grid%truelat1, grid%truelat2, &
276 coarse_ix, coarse_jy, start_x, start_y)
277 end if
278 end if
279 ! [8.8] Write Ascii radiance OMB and OMA file
280
281 if (write_oa_rad_ascii) then
282 write(unit=stdout,fmt='(A)') 'Writing radiance OMB and OMA ascii file'
283 write(unit=stdout,fmt=*) " "
284 call da_write_oa_rad_ascii (grid%xp,ob,iv,re)
285 end if
286
287 !------------------------------------------------------------------------
288 ! [8.0] Output WRFVAR analysis and analysis increments:
289 !------------------------------------------------------------------------
290
291 call da_transfer_xatoanalysis( it, xbx, grid, config_flags)
292 end do
293
294 !---------------------------------------------------------------------------
295 ! [9.0] Tidy up:
296 !---------------------------------------------------------------------------
297
298 if (var4d) then
299 call da_system("touch wrf_stop_now")
300 end if
301
302 deallocate (cvt)
303 deallocate (xhat)
304
305 if (use_rad) then
306 do i =1, iv%num_inst
307 deallocate (j % jo % rad(i) % jo_ichan)
308 deallocate (j % jo % rad(i) % num_ichan)
309 deallocate (satinfo(i) % ichan)
310 deallocate (satinfo(i) % iuse)
311 deallocate (satinfo(i) % error)
312 deallocate (satinfo(i) % polar)
313
314 if (read_biascoef) then
315 deallocate (satinfo(i) % scanbias)
316 deallocate (satinfo(i) % scanbias_b)
317 deallocate (satinfo(i) % bcoef)
318 deallocate (satinfo(i) % bcoef0)
319 deallocate (satinfo(i) % error_std)
320 end if
321
322 deallocate (ob%instid(i) % ichan)
323 deallocate (iv%instid(i)%ichan)
324
325 if (iv%instid(i)%num_rad > 0) then
326 deallocate (iv%instid(i)%info)
327 deallocate (iv%instid(i)%loc)
328 deallocate (iv%instid(i)%loc_i)
329 deallocate (iv%instid(i)%loc_j)
330 deallocate (iv%instid(i)%loc_k)
331 deallocate (iv%instid(i)%loc_dx)
332 deallocate (iv%instid(i)%loc_dy)
333 deallocate (iv%instid(i)%loc_dz)
334 deallocate (iv%instid(i)%loc_dxm)
335 deallocate (iv%instid(i)%loc_dym)
336 deallocate (iv%instid(i)%loc_dzm)
337 deallocate (iv%instid(i)%zk)
338 deallocate (iv%instid(i)%t)
339 deallocate (iv%instid(i)%mr)
340 deallocate (iv%instid(i)%tm)
341 deallocate (iv%instid(i)%qm)
342 deallocate (iv%instid(i)%qrn)
343 deallocate (iv%instid(i)%qcw)
344 deallocate (iv%instid(i)%qci)
345 deallocate (iv%instid(i)%qsn)
346 deallocate (iv%instid(i)%qgr)
347 deallocate (iv%instid(i)%pm)
348 deallocate (iv%instid(i)%pf)
349 deallocate (iv%instid(i)%u10)
350 deallocate (iv%instid(i)%v10)
351 deallocate (iv%instid(i)%t2m)
352 deallocate (iv%instid(i)%q2m)
353 deallocate (iv%instid(i)%mr2m)
354 deallocate (iv%instid(i)%psfc)
355 deallocate (iv%instid(i)%ts)
356 deallocate (iv%instid(i)%smois)
357 deallocate (iv%instid(i)%tslb)
358 deallocate (iv%instid(i)%snowh)
359 deallocate (iv%instid(i)%isflg)
360 deallocate (iv%instid(i)%soiltyp)
361 deallocate (iv%instid(i)%landsea_mask)
362 deallocate (iv%instid(i)%elevation)
363 deallocate (iv%instid(i)%vegfra)
364 deallocate (iv%instid(i)%vegtyp)
365 deallocate (iv%instid(i)%clwp)
366 deallocate (iv%instid(i)%ps)
367 deallocate (iv%instid(i)%tb_xb)
368 deallocate (iv%instid(i)%tb_qc)
369 deallocate (iv%instid(i)%tb_inv)
370 deallocate (iv%instid(i)%tb_error)
371 deallocate (iv%instid(i)%emiss)
372 deallocate (iv%instid(i)%scanpos)
373 deallocate (iv%instid(i)%scanline)
374 deallocate (iv%instid(i)%ifgat)
375 deallocate (iv%instid(i)%cloud_flag)
376 deallocate (iv%instid(i)%satzen)
377 deallocate (iv%instid(i)%satazi)
378 deallocate (iv%instid(i)%solzen)
379 deallocate (iv%instid(i)%solazi)
380 deallocate (iv%instid(i)%proc_domain)
381
382 if (rtm_option == rtm_option_crtm) then
383 deallocate(iv%instid(i)%water_coverage)
384 deallocate(iv%instid(i)%land_coverage)
385 deallocate(iv%instid(i)%ice_coverage)
386 deallocate(iv%instid(i)%snow_coverage)
387 !deallocate(iv%instid(i)%ps_jacobian)
388 !deallocate(iv%instid(i)%t_jacobian)
389 !deallocate(iv%instid(i)%q_jacobian)
390 end if
391 end if
392 #ifdef RTTOV
393 if (rtm_option == rtm_option_rttov) then
394 call rttov_dealloc_coef (ierr,coefs(i))
395 end if
396 #endif
397 end do
398 deallocate (iv%instid)
399 deallocate (j % jo % rad)
400 deallocate (satinfo)
401 deallocate (time_slots)
402 #ifdef RTTOV
403 if (rtm_option == rtm_option_rttov) then
404 deallocate (coefs)
405 end if
406 #endif
407 end if
408
409 call da_deallocate_observations(iv)
410 call da_deallocate_y (re)
411 call da_deallocate_y (y)
412 call da_deallocate_y (ob)
413 call da_deallocate_background_errors (be)
414
415 if (xbx%pad_num > 0) then
416 deallocate(xbx%pad_loc)
417 deallocate(xbx%pad_pos)
418 end if
419
420 deallocate (xbx % fft_factors_x)
421 deallocate (xbx % fft_factors_y)
422 deallocate (xbx % fft_coeffs)
423 deallocate (xbx % trig_functs_x)
424 deallocate (xbx % trig_functs_y)
425
426 if (global) then
427 deallocate (xbx%coslat)
428 deallocate (xbx%sinlat)
429 deallocate (xbx%coslon)
430 deallocate (xbx%sinlon)
431 deallocate (xbx%int_wgts)
432 deallocate (xbx%alp)
433 deallocate (xbx%wsave)
434 if (grid%xb%jts == grid%xb%jds) then
435 deallocate(cos_xls)
436 deallocate(sin_xls)
437 end if
438
439 if (grid%xb%jte == grid%xb%jde) then
440 deallocate(cos_xle)
441 deallocate(sin_xle)
442 end if
443 end if
444
445 deallocate (xbx % latc_mean)
446
447 #ifdef DM_PARALLEL
448 call mpi_barrier(comm,ierr)
449 #endif
450
451 if (trace_use) call da_trace_exit("da_solve")
452
453 contains
454
455 #include "da_solve_init.inc"
456
457 end subroutine da_solve
458