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) 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 end if
167
168 ! [8.2] Calculate innovation vector (O-B):
169
170 call da_get_innov_vector (it, ob, iv, grid , config_flags)
171
172 if (test_transforms .or. test_wrfvar) then
173 call da_check (grid, config_flags, cv_size, xbx, be, grid%ep, iv, &
174 grid%vv, grid%vp, y)
175 if (var4d) then
176 call da_system ("touch wrf_stop_now")
177 end if
178 ! No point continuing, as data corrupted
179 call wrf_shutdown
180 stop
181 end if
182
183 ! [8.4] Minimize cost function:
184
185 call da_allocate_y (iv, re)
186 call da_allocate_y (iv, y)
187
188 call da_minimise_cg( grid, config_flags, it, be % cv % size, &
189 xbx, be, iv, j_grad_norm_target, xhat, cvt, re, y, j)
190
191 !------------------------------------------------------------------------
192
193 ! [8.5] Update latest analysis solution:
194
195 call da_transform_vtox (grid,cv_size,xbx,be,grid%ep,xhat,grid%vv,grid%vp)
196
197 ! [8.6] Only when use_radarobs = .false. and calc_w_increment =.true.,
198 ! the w_increment need to be diagnosed:
199
200 if (calc_w_increment .and. .not. use_radarobs) then
201 call da_uvprho_to_w_lin (grid)
202
203 #ifdef DM_PARALLEL
204 #include "HALO_RADAR_XA_W.inc"
205 #endif
206 end if
207
208 ! [8.7] Write out diagnostics
209
210 call da_write_diagnostics (grid, ob, iv, re, y, j)
211
212 ! Write "clean" QCed observations if requested:
213 if (anal_type_qcobs) then
214 if (it == 1) then
215 call da_write_filtered_obs (grid, ob, iv, &
216 coarse_ix, coarse_jy, start_x, start_y)
217 end if
218 end if
219
220 ! [8.3] Interpolate x_g to low resolution grid
221
222 ! [8.8] Write Ascii radiance OMB and OMA file
223
224 if (write_oa_rad_ascii) then
225 write(unit=stdout,fmt='(A)') 'Writing radiance OMB and OMA ascii file'
226 write(unit=stdout,fmt=*) " "
227 call da_write_oa_rad_ascii (ob,iv,re)
228 end if
229
230 !------------------------------------------------------------------------
231 ! [8.0] Output WRFVAR analysis and analysis increments:
232 !------------------------------------------------------------------------
233
234 call da_transfer_xatoanalysis (it, xbx, grid, config_flags)
235 end do
236
237 !---------------------------------------------------------------------------
238 ! [9.0] Tidy up:
239 !---------------------------------------------------------------------------
240
241 if (var4d) then
242 call da_system ("touch wrf_stop_now")
243 end if
244
245 deallocate (cvt)
246 deallocate (xhat)
247
248 if (use_rad) then
249 do i =1, iv%num_inst
250 deallocate (j % jo % rad(i) % jo_ichan)
251 deallocate (j % jo % rad(i) % num_ichan)
252 deallocate (satinfo(i) % ichan)
253 deallocate (satinfo(i) % iuse)
254 deallocate (satinfo(i) % error)
255 deallocate (satinfo(i) % polar)
256
257 if (read_biascoef) then
258 deallocate (satinfo(i) % scanbias)
259 deallocate (satinfo(i) % scanbias_b)
260 deallocate (satinfo(i) % bcoef)
261 deallocate (satinfo(i) % bcoef0)
262 deallocate (satinfo(i) % error_std)
263 end if
264
265 deallocate (ob%instid(i) % ichan)
266 deallocate (iv%instid(i)%ichan)
267
268 if (iv%instid(i)%num_rad > 0) then
269
270 deallocate (iv%instid(i)%info%levels)
271 deallocate (iv%instid(i)%info%lat)
272 deallocate (iv%instid(i)%info%lon)
273 deallocate (iv%instid(i)%info%elv)
274
275 deallocate (iv%instid(i)%info%i)
276 deallocate (iv%instid(i)%info%j)
277 deallocate (iv%instid(i)%info%k)
278 deallocate (iv%instid(i)%info%zk)
279 deallocate (iv%instid(i)%info%dx)
280 deallocate (iv%instid(i)%info%dy)
281 deallocate (iv%instid(i)%info%dz)
282 deallocate (iv%instid(i)%info%dxm)
283 deallocate (iv%instid(i)%info%dym)
284 deallocate (iv%instid(i)%info%dzm)
285
286 deallocate (iv%instid(i)%t)
287 deallocate (iv%instid(i)%mr)
288 deallocate (iv%instid(i)%tm)
289 deallocate (iv%instid(i)%qm)
290 deallocate (iv%instid(i)%qrn)
291 deallocate (iv%instid(i)%qcw)
292 deallocate (iv%instid(i)%qci)
293 deallocate (iv%instid(i)%qsn)
294 deallocate (iv%instid(i)%qgr)
295 deallocate (iv%instid(i)%pm)
296 deallocate (iv%instid(i)%pf)
297 deallocate (iv%instid(i)%u10)
298 deallocate (iv%instid(i)%v10)
299 deallocate (iv%instid(i)%t2m)
300 deallocate (iv%instid(i)%q2m)
301 deallocate (iv%instid(i)%mr2m)
302 deallocate (iv%instid(i)%psfc)
303 deallocate (iv%instid(i)%ts)
304 deallocate (iv%instid(i)%smois)
305 deallocate (iv%instid(i)%tslb)
306 deallocate (iv%instid(i)%snowh)
307 deallocate (iv%instid(i)%isflg)
308 deallocate (iv%instid(i)%soiltyp)
309 deallocate (iv%instid(i)%landsea_mask)
310 deallocate (iv%instid(i)%elevation)
311 deallocate (iv%instid(i)%vegfra)
312 deallocate (iv%instid(i)%vegtyp)
313 deallocate (iv%instid(i)%clwp)
314 deallocate (iv%instid(i)%ps)
315 deallocate (iv%instid(i)%tb_xb)
316 deallocate (iv%instid(i)%tb_qc)
317 deallocate (iv%instid(i)%tb_inv)
318 deallocate (iv%instid(i)%tb_error)
319 deallocate (iv%instid(i)%emiss)
320 deallocate (iv%instid(i)%scanpos)
321 deallocate (iv%instid(i)%scanline)
322 deallocate (iv%instid(i)%ifgat)
323 deallocate (iv%instid(i)%cloud_flag)
324 deallocate (iv%instid(i)%satzen)
325 deallocate (iv%instid(i)%satazi)
326 deallocate (iv%instid(i)%solzen)
327 deallocate (iv%instid(i)%solazi)
328
329 if (rtm_option == rtm_option_crtm) then
330 deallocate(iv%instid(i)%water_coverage)
331 deallocate(iv%instid(i)%land_coverage)
332 deallocate(iv%instid(i)%ice_coverage)
333 deallocate(iv%instid(i)%snow_coverage)
334 ! deallocate(iv%instid(i)%ps_jacobian)
335 ! deallocate(iv%instid(i)%t_jacobian)
336 ! deallocate(iv%instid(i)%q_jacobian)
337 end if
338 end if
339 #ifdef RTTOV
340 if (rtm_option == rtm_option_rttov) then
341 call rttov_dealloc_coef (ierr,coefs(i))
342 end if
343 #endif
344 end do
345 deallocate (iv%instid)
346 deallocate (j % jo % rad)
347 deallocate (satinfo)
348 deallocate (time_slots)
349 #ifdef RTTOV
350 if (rtm_option == rtm_option_rttov) then
351 deallocate (coefs)
352 end if
353 #endif
354 end if
355
356 call da_deallocate_observations (iv)
357 call da_deallocate_y (re)
358 call da_deallocate_y (y)
359 call da_deallocate_y (ob)
360 call da_deallocate_background_errors (be)
361
362 if (xbx%pad_num > 0) then
363 deallocate (xbx%pad_loc)
364 deallocate (xbx%pad_pos)
365 end if
366
367 deallocate (xbx % fft_factors_x)
368 deallocate (xbx % fft_factors_y)
369 deallocate (xbx % fft_coeffs)
370 deallocate (xbx % trig_functs_x)
371 deallocate (xbx % trig_functs_y)
372
373 if (global) then
374 deallocate (xbx%coslat)
375 deallocate (xbx%sinlat)
376 deallocate (xbx%coslon)
377 deallocate (xbx%sinlon)
378 deallocate (xbx%int_wgts)
379 deallocate (xbx%alp)
380 deallocate (xbx%wsave)
381 if (jts == jds) then
382 deallocate (cos_xls)
383 deallocate (sin_xls)
384 end if
385
386 if (jte == jde) then
387 deallocate (cos_xle)
388 deallocate (sin_xle)
389 end if
390 end if
391
392 deallocate (xbx % latc_mean)
393
394 #ifdef DM_PARALLEL
395 call mpi_barrier (comm,ierr)
396 #endif
397
398 if (trace_use) call da_trace_exit ("da_solve")
399
400 contains
401
402 #include "da_solve_init.inc"
403
404 end subroutine da_solve
405