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