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