var/da/da_main/da_solve.inc

   1   subroutine da_solve ( grid , config_flags)
   2
   3    !-----------------------------------------------------------------------
   4    ! Purpose: TBD
   5    !
   6    ! Edited 09/06/2012: Allow for variable ntmax for each outer loop (Mike Kavulich)
   7    !-----------------------------------------------------------------------
   8
   9    implicit none
  10
  11
  12
  13    type (domain),               intent(inout) :: grid
  14    type (grid_config_rec_type), intent(inout) :: config_flags
  15
  16    type (xbx_type)              :: xbx         ! For header & non-grid arrays.
  17    type (be_type)               :: be          ! Background error structure.
  18    real, allocatable            :: cvt(:)      ! Control variable structure.
  19    real, allocatable            :: xhat(:)     ! Control variable structure.
  20    real, allocatable            :: qhat(:,:)   ! Control variable structure.
  21    real*8, allocatable          :: eignvec(:,:)
  22    real*8, allocatable          :: eignval(:)
  23 !   real, allocatable            :: full_eignvec(:)
  24    type (y_type)                :: ob          ! Observation structure.
  25    type (iv_type)               :: iv          ! Obs. increment structure.
  26    type (y_type)                :: re          ! Residual (o-a) structure.
  27    type (y_type)                :: y           ! y = H(x_inc) structure.
  28    integer                      :: it          ! External loop counter.
  29    integer                      :: neign
  30    type (j_type)                :: j           ! Cost function.
  31    type (y_type)                :: jo_grad_y   ! Grad_y(jo)
  32
  33    integer                      :: cv_size, i, ichan, k, n
  34    real                         :: j_grad_norm_target ! Target j norm.
  35
  36    character (len=3)            :: ci
  37    character (len=2)            :: outerloop
  38    character (len=256)          :: timestr
  39    integer   :: min_yyyy,min_mm,min_dd,min_hh,min_mn,min_ss
  40    integer   :: max_yyyy,max_mm,max_dd,max_hh,max_mn,max_ss
  41    character :: s
  42    real*8    :: time_min, time_max
  43    integer   :: jl_start, jl_end
  44    character(len=256) :: timestr1
  45    type(x_type) :: shuffle
  46
  47    real, allocatable :: grid_box_area(:,:), mapfac(:,:)
  48    real, allocatable :: v1(:,:,:),v2(:,:,:),v3(:,:,:),v4(:,:,:),v5(:,:,:)
  49    real, allocatable :: v6(:,:,:),v7(:,:,:),v8(:,:,:),v9(:,:,:),v10(:,:,:),v11(:,:,:)
  50    character (len=10)    :: variable_name
  51
  52    integer :: iwin, num_subtwindow
  53    integer :: ni, nj, ii, jj
  54    real    :: dx, dy, dxm, dym, cxp0, cxp1, cyp0, cyp1
  55    real,    external :: nest_loc_of_cg  ! from share/interp_fcn.F
  56    integer, external :: compute_CGLL    ! from share/interp_fcn.F
  57
  58    integer   :: vp_unit, iost
  59    character(len=13) :: vpfile ! vp_input.0001
  60    integer   :: i1,i2,i3,i4,i5,i6
  61    logical :: ex
  62    logical :: offline_varbc
  63
  64    character(len=10) :: this_time
  65    character(len=8)  :: this_date
  66
  67    if (trace_use) call da_trace_entry("da_solve")
  68
  69 #ifdef DM_PARALLEL
  70    call mpi_barrier(comm,ierr)
  71 #endif
  72
  73    if ( config_flags%use_baseparam_fr_nml ) then
  74       call nl_get_base_pres  ( 1 , base_pres )
  75       call nl_get_base_temp  ( 1 , base_temp )
  76       call nl_get_base_lapse ( 1 , base_lapse )
  77       call nl_get_iso_temp   ( 1 , iso_temp )
  78       if ( iso_temp .NE. grid%tiso ) THEN
  79          write(unit=message(1),fmt='(A)') &
  80            'Namelist iso_temp does not equal iso_temp from fg. Reset namelist value and rerun.'
  81          write(unit=message(2),fmt='(A,F10.5)')'Namelist iso_temp   = ',iso_temp
  82          write(unit=message(3),fmt='(A,F10.5)')'Background iso_temp = ',grid%tiso
  83          call da_error(__FILE__,__LINE__,message(1:3))
  84       end if
  85       call nl_get_base_pres_strat  ( 1, base_pres_strat )
  86       call nl_get_base_lapse_strat ( 1, base_lapse_strat )
  87
  88       grid%p00   = base_pres
  89       grid%t00   = base_temp
  90       grid%tlp   = base_lapse
  91       grid%tiso  = iso_temp
  92       grid%p_strat   = base_pres_strat
  93       grid%tlp_strat = base_lapse_strat
  94    else
  95       base_pres  = grid%p00
  96       base_temp  = grid%t00
  97       base_lapse = grid%tlp
  98       iso_temp   = grid%tiso
  99       base_pres_strat  = grid%p_strat
 100       base_lapse_strat = grid%tlp_strat
 101       if ( base_temp < 100.0 .or. base_pres < 10000.0 ) then
 102          write(unit=message(1),fmt='(A)') &
 103          'did not find base state parameters in fg. Add use_baseparam_fr_nml = .t. in &dynamics and rerun '
 104          call da_error(__FILE__,__LINE__,message(1:1))
 105       end if
 106    end if
 107
 108    ! Calculate the num_fgat_time based on time_window_min, time_window_max
 109     if ( var4d ) then
 110        if (time_step == 0) then
 111           write(unit=message(1),fmt='(A)') &
 112           'For 4DVAR, in the &domains namelist, "time_step" must be set to a non-zero value'
 113           call da_error(__FILE__,__LINE__,message(1:1))
 114        endif
 115        read(unit=time_window_min,fmt='(i4,5(a1,i2))') min_yyyy,s,min_mm,s,min_dd,s,min_hh,s,min_mn,s,min_ss
 116        read(unit=time_window_max,fmt='(i4,5(a1,i2))') max_yyyy,s,max_mm,s,max_dd,s,max_hh,s,max_mn,s,max_ss
 117        call da_get_julian_time(min_yyyy,min_mm,min_dd,min_hh,min_mn,time_min)
 118        call da_get_julian_time(max_yyyy,max_mm,max_dd,max_hh,max_mn,time_max)
 119        if ( var4d_bin < time_step ) call nl_set_var4d_bin (1, time_step)
 120        time_max = (time_max - time_min) * 60   ! unit is : seconds
 121        num_fgat_time = NINT(time_max/var4d_bin)
 122        if ( NINT(time_max/var4d_bin)*var4d_bin .ne. NINT(time_max) ) then
 123           write(unit=message(1),fmt='(A)') &
 124           '4DVAR assimilation window must be evenly divisible by var4d_bin!'
 125           write(unit=message(2),fmt='(A,I7)') &
 126           'var4d_bin       = ',var4d_bin
 127           write(unit=message(3),fmt='(A,A)') &
 128           'time_window_max = ',time_window_max
 129           write(unit=message(4),fmt='(A,A)') &
 130           'time_window_min = ',time_window_min
 131           write(unit=message(5),fmt='(A,F10.4)') &
 132           'time_window_max - time_window_min = ',time_max
 133           write(unit=message(6),fmt='(A)')'Change var4d_bin, time_window_max, or time_window_min in namelist and rerun'
 134           call da_error(__FILE__,__LINE__,message(1:6))
 135        endif
 136        if ( var4d_bin/time_step*time_step .ne. var4d_bin ) then
 137           write(unit=message(1),fmt='(A)') &
 138           'var4d_bin must be evenly divisible by time_step!'
 139           write(unit=message(2),fmt='(A,I7)') &
 140           'var4d_bin = ',var4d_bin
 141           write(unit=message(3),fmt='(A,I7)') &
 142           'time_step = ',time_step
 143           write(unit=message(4),fmt='(A)')'Change var4d_bin or time_step in namelist and rerun'
 144           call da_error(__FILE__,__LINE__,message(1:4))
 145        endif
 146
 147        num_fgat_time = num_fgat_time + 1
 148        write(unit=message(1),fmt='(a,i10)') 'num_fgat_time is: ', num_fgat_time
 149        call da_message(message(1:1))
 150        if ( use_rainobs ) then
 151           allocate (fgat_rain_flags(1:num_fgat_time))
 152           fgat_rain_flags = .false.
 153           if ( INT(var4d_bin_rain/var4d_bin)*var4d_bin .ne. INT(var4d_bin_rain) ) then
 154              write(unit=message(1),fmt='(A,A,2I7)') &
 155              'Please change var4d_bin_rain in namelist and rerun==>', 'var4d_bin_rain, var4d_bin:',var4d_bin_rain,var4d_bin
 156              call da_error(__FILE__,__LINE__,message(1:1))
 157           endif
 158           do n = 1, num_fgat_time, INT(var4d_bin_rain/var4d_bin)
 159              fgat_rain_flags(n) = .true.
 160           end do
 161        end if
 162     endif
 163
 164    !---------------------------------------------------------------------------
 165    ! [1.0] Initial checks
 166    !---------------------------------------------------------------------------
 167
 168    if (cv_options_hum /= cv_options_hum_specific_humidity .and. &
 169        cv_options_hum /= cv_options_hum_relative_humidity) then
 170       write(unit=message(1),fmt='(A,I3)') &
 171          'Invalid cv_options_hum = ', cv_options_hum
 172       call da_error(__FILE__,__LINE__,message(1:1))
 173    end if
 174
 175    if (vert_corr == vert_corr_2) then
 176       if (vertical_ip < vertical_ip_0 .or. vertical_ip > vertical_ip_delta_p) then
 177          write (unit=message(1),fmt='(A,I3)') &
 178            'Invalid vertical_ip = ', vertical_ip
 179          call da_error(__FILE__,__LINE__,message(1:1))
 180       end if
 181    end if
 182
 183    if( use_rf )then
 184       if (0.5 * real(rf_passes) /= real(rf_passes / 2)) then
 185          write(unit=stdout,fmt='(A,I4,A)') &
 186             'rf_passes = ', rf_passes, ' .Should be even.'
 187          rf_passes = int(real(rf_passes / 2))
 188          write(unit=stdout,fmt='(A,I4)') 'Resetting rf_passes = ', rf_passes
 189       end if
 190    else
 191       write(stdout,'("da_solve: using wavelet transform")')
 192    endif
 193
 194    if ( anal_type_hybrid_dual_res .and. alphacv_method .ne. alphacv_method_xa ) then
 195       write (unit=message(1),fmt='(A)') &
 196         'Dual-res hybrid only with alphacv_method = 2'
 197       call da_error(__FILE__,__LINE__,message(1:1))
 198    endif
 199
 200    if (anal_type_randomcv) then
 201       write_gts_omb_oma  = .false.
 202       write_rej_obs_conv = .false.
 203       write_unpert_obs   = .false.
 204    end if
 205
 206    if ( cv_options == 3 ) then
 207       if ( cloud_cv_options > 0 ) then
 208          write(unit=message(1),fmt='(A)') &
 209             'Cloud control variables are not implemented for cv_options = 3'
 210          write(unit=message(2),fmt='(A)') &
 211             'Resetting cloud_cv_options = 0 for cv_options = 3'
 212          cloud_cv_options = 0
 213          call da_warning(__FILE__,__LINE__,message(1:2))
 214       end if
 215       if ( ensdim_alpha > 0 ) then
 216          write(unit=message(1),fmt='(A)') &
 217             'Alpha control variables are not implemented for cv_options = 3'
 218          write(unit=message(2),fmt='(A)') &
 219             'Please set cv_options=5 or 7, or ensdim_alpha=0'
 220          call da_error(__FILE__,__LINE__,message(1:2))
 221       end if
 222    end if
 223
 224    if ( use_radarobs ) then
 225       if ( use_radar_rf .and. use_radar_rhv ) then
 226          write(unit=message(1),fmt='(A)') "Both 'use_radar_rf' and 'use_radar_rhv' are set to true"
 227          write(unit=message(2),fmt='(A)') "You must only choose one of these options"
 228          call da_error(__FILE__,__LINE__,message(1:2))
 229       end if
 230       if ( (use_radar_rf .and. radar_rf_opt==1) .and. cloud_cv_options /= 1 ) then
 231          write(unit=message(1),fmt='(A)') "Please set cloud_cv_options=1 for use_radar_rf and radar_rf_opt=1"
 232          call da_error(__FILE__,__LINE__,message(1:1))
 233       end if
 234       if ( use_radar_rhv .and. cloud_cv_options == 1 ) then
 235          write(unit=message(1),fmt='(A)') "Please set cloud_cv_options=3 or 2 (2 requires cloudy be.dat) for use_radar_rhv"
 236          call da_error(__FILE__,__LINE__,message(1:1))
 237       end if
 238    end if
 239
 240    if ( ensdim_alpha > 0 .and. alpha_hydrometeors ) then
 241       if ( cloud_cv_options == 1 ) then
 242          write(unit=message(1),fmt='(A)') "alpha_hydrometeors is not implemented for cloud_cv_options=1"
 243          write(unit=message(2),fmt='(A)') "Please choose cloud_cv_options=3"
 244          call da_error(__FILE__,__LINE__,message(1:2))
 245       end if
 246    end if
 247
 248    if ( use_gpsrefobs .and. use_gpsephobs ) then
 249       write(unit=message(1),fmt='(A,I3)') &
 250          'Both use_gpsrefobs and use_gpsephobs are true. Should choose either use_gpsrefobs or use_gpsephobs'
 251       call da_error(__FILE__,__LINE__,message(1:1))
 252    end if
 253
 254    !---------------------------------------------------------------------------
 255    ! [2.0] Initialise wrfvar parameters:
 256    !---------------------------------------------------------------------------
 257
 258    if ( anal_type_hybrid_dual_res ) then
 259
 260       !---------------------------------
 261       ! Get full ensemble grid dimensions
 262       !---------------------------------
 263       call da_solve_init(ensemble_grid &
 264 #include "actual_new_args.inc"
 265 )
 266       ide_ens = ide ! these are unstaggered dimensions of the full ensemble domain
 267       jde_ens = jde
 268       kde_ens = kde
 269
 270       !---------------------------------------
 271       ! Get "intermediate" grid sizes and tiles
 272       !---------------------------------------
 273       call da_solve_init(grid%intermediate_grid &
 274 #include "actual_new_args.inc"
 275 )
 276
 277       ! these are unstaggered dimensions of the "intermediate" ensemble domain
 278       !  The intermediate grid is the coarse (ensemble) domain that is co-located with the
 279       !  hi-resolution (analysis) grid
 280
 281       ids_int = ids ; jds_int = jds ; kds_int = kds
 282       ide_int = ide ; jde_int = jde ; kde_int = kde
 283
 284       its_int = its ; ite_int = ite
 285       jts_int = jts ; jte_int = jte
 286       kts_int = kts ; kte_int = kte
 287
 288       ims_int = ims ; ime_int = ime
 289       jms_int = jms ; jme_int = jme
 290       kms_int = kms ; kme_int = kme
 291
 292       ips_int = ips ; ipe_int = ipe
 293       jps_int = jps ; jpe_int = jpe
 294       kps_int = kps ; kpe_int = kpe
 295
 296
 297       grid%imask_xstag = 1   ; grid%imask_ystag = 1
 298       grid%imask_nostag = 1  ; grid%imask_xystag = 1
 299
 300       !---------------------------------------------------------------------------
 301       ! De-allocate parts of grid and replace with grid%intermediate_grid dimensions
 302       !---------------------------------------------------------------------------
 303       call reallocate_analysis_grid(grid)
 304
 305       !----------------------------------------------------------
 306       ! Allocate and initialize some of grid%intermediate_grid
 307       !----------------------------------------------------------
 308       call allocate_intermediate_grid(grid%intermediate_grid)
 309
 310
 311       !---------------------------------------
 312       ! Get map projection information for the ensemble
 313       !---------------------------------------
 314
 315        call da_setup_firstguess(xbx, ensemble_grid, config_flags, .true. )
 316        map_info_ens = map_info ! map_info is read in from da_tools.f90, call it something else
 317
 318    endif
 319
 320    call da_solve_init(grid &
 321 #include "actual_new_args.inc"
 322 )
 323
 324    if ( .not. anal_type_hybrid_dual_res ) then
 325       ide_ens = ide ; jde_ens = jde ; kde_ens = kde
 326
 327       ids_int = ids ; ide_int = ide
 328       jds_int = jds ; jde_int = jde
 329       kds_int = kds ; kde_int = kde
 330
 331       its_int = its ; ite_int = ite
 332       jts_int = jts ; jte_int = jte
 333       kts_int = kts ; kte_int = kte
 334
 335       ims_int = ims ; ime_int = ime
 336       jms_int = jms ; jme_int = jme
 337       kms_int = kms ; kme_int = kme
 338
 339       ips_int = ips ; ipe_int = ipe
 340       jps_int = jps ; jpe_int = jpe
 341       kps_int = kps ; kpe_int = kpe
 342    endif
 343
 344    !---------------------------------------------------------------------------
 345    ! [3.0] Set up first guess field (grid%xb):
 346    !---------------------------------------------------------------------------
 347
 348    ! for WRF hybrid coordinate
 349    allocate ( c1f(kms:kme) )
 350    allocate ( c2f(kms:kme) )
 351    allocate ( c3f(kms:kme) )
 352    allocate ( c4f(kms:kme) )
 353    allocate ( c1h(kms:kme) )
 354    allocate ( c2h(kms:kme) )
 355    allocate ( c3h(kms:kme) )
 356    allocate ( c4h(kms:kme) )
 357    c1f = grid%c1f
 358    c2f = grid%c2f
 359    c3f = grid%c3f
 360    c4f = grid%c4f
 361    c1h = grid%c1h
 362    c2h = grid%c2h
 363    c3h = grid%c3h
 364    c4h = grid%c4h
 365    ! For WRFDA, namelist hybrid_opt setting is ignored.
 366    ! grid%hybrid_opt is set in da_med_initialdata_input.inc based on
 367    ! HYBRID_OPT global attribute contained in the input file.
 368    ! For input file prior to V3.9, grid%hybrid_opt is set to 0.
 369    if ( grid%hybrid_opt <= 0 ) then
 370       c3f = grid%znw
 371       c3h = grid%znu
 372       c4f = 0.0
 373       c4h = 0.0
 374       c1f = 1.0
 375       c1h = 1.0
 376       c2f = 0.0
 377       c2h = 0.0
 378    end if
 379
 380    call da_setup_firstguess(xbx, grid, config_flags, .false.)
 381
 382    if ( anal_type_hybrid_dual_res ) then
 383
 384       !
 385       ! Get ensemble grid mapfactor on entire coarse grid
 386       !
 387       variable_name = 'MAPFAC_M'
 388
 389       allocate( grid_box_area(1:ide_ens,1:jde_ens), mapfac(1:ide_ens,1:jde_ens) )
 390       call da_get_var_2d_real_cdf( input_file_ens, variable_name, mapfac, ide_ens, jde_ens, 1, .false. )
 391       grid_box_area(:,:) = ( (ensemble_grid%dx)/mapfac(:,:) )**2
 392       grid%intermediate_grid%xb%grid_box_area(its_int:ite_int,jts_int:jte_int) = grid_box_area(its_int:ite_int,jts_int:jte_int)
 393       deallocate(mapfac,grid_box_area)
 394
 395       !
 396       ! get fine grid locs relative to coarse grid for bilinear interplation
 397       ! between coarse (ensemble) and fine (analysis) grids.
 398       ! use functions compute_CGLL and nest_loc_of_cg taken from
 399       ! share/interp_fcn.F
 400       !
 401       allocate(aens_locs(its:ite, jts:jte)) ! From da_control
 402       do nj = jts, jte
 403          jj = compute_CGLL (nj, grid%j_parent_start, grid%parent_grid_ratio, 1)
 404          cyp0 = nest_loc_of_cg (jj  , grid%j_parent_start, grid%parent_grid_ratio, 0)
 405          cyp1 = nest_loc_of_cg (jj+1, grid%j_parent_start, grid%parent_grid_ratio, 0)
 406          dym = ( cyp1 - real(nj) ) / ( cyp1 - cyp0 )
 407          dy = 1.0 - dym
 408          do ni = its, ite
 409             ii = compute_CGLL (ni, grid%i_parent_start, grid%parent_grid_ratio, 1)
 410             cxp0 = nest_loc_of_cg (ii  , grid%i_parent_start, grid%parent_grid_ratio, 0)
 411             cxp1 = nest_loc_of_cg (ii+1, grid%i_parent_start, grid%parent_grid_ratio, 0)
 412             dxm = ( cxp1 - real(ni) ) / ( cxp1 - cxp0 )
 413             dx = 1.0 - dxm
 414             aens_locs(ni,nj)%i   = ii
 415             aens_locs(ni,nj)%j   = jj
 416             aens_locs(ni,nj)%dx  = dx
 417             aens_locs(ni,nj)%dy  = dy
 418             aens_locs(ni,nj)%dxm = dxm
 419             aens_locs(ni,nj)%dym = dym
 420          end do
 421       end do
 422
 423    end if !anal_type_hybrid_dual_res
 424
 425    !---------------------------------------------------------------------------
 426    ! [4.0] Set up observations (ob):
 427    !---------------------------------------------------------------------------
 428    call da_setup_obs_structures (grid, ob, iv, j)
 429    if (use_rad) then
 430       allocate (j % jo % rad(1:iv%num_inst))
 431       do i=1,iv%num_inst
 432          allocate (j % jo % rad(i) % jo_ichan(iv%instid(i)%nchan))
 433          allocate (j % jo % rad(i) % num_ichan(iv%instid(i)%nchan))
 434       end do
 435    end if
 436
 437    !---------------------------------------------------------------------------
 438    ! [4.1] Observer (ANAL_TYPE="VERIFY")
 439    !---------------------------------------------------------------------------
 440
 441    if (anal_type_verify) then
 442       check_max_iv = .false.
 443       ntmax=0
 444       it = 1
 445       num_qcstat_conv=0
 446
 447 #if defined(RTTOV) || defined(CRTM)
 448       if (use_rad .and. (use_varbc.or.freeze_varbc)) call da_varbc_init(iv, be)
 449 #endif
 450
 451       call da_get_innov_vector (it, num_qcstat_conv, ob, iv, grid , config_flags)
 452       call da_allocate_y (iv, re)
 453
 454       ! write out O-B statistics
 455       call da_write_diagnostics(it, grid,num_qcstat_conv, ob, iv, re, y, j)
 456
 457       ! write out Gradient of Jo for adjoint sensitivity
 458       if (adj_sens) then
 459          cv_size = 1
 460          allocate (xhat(cv_size))
 461          call da_allocate_y (iv, y)
 462          call da_allocate_y (iv, jo_grad_y)
 463
 464          call da_calculate_residual(iv, y, re)
 465          call da_calculate_grady(iv, re, jo_grad_y)
 466          call da_zero_x(grid%xa)
 467
 468          call da_transform_xtoy_adj(cv_size, xhat, grid, iv, jo_grad_y, grid%xa)
 469          call da_transform_xtoxa_adj(grid)
 470          call da_transfer_wrftltoxa_adj(grid, config_flags, 'fcst', timestr)
 471
 472          call da_deallocate_y (y)
 473          call da_deallocate_y (jo_grad_y)
 474       end if
 475
 476       call da_deallocate_y(re)
 477       call da_deallocate_observations (iv)
 478       if (trace_use) call da_trace_exit ("da_solve")
 479       return
 480    end if
 481
 482    !---------------------------------------------------------------------------
 483    ! [5.0] Set up control variable:
 484    !---------------------------------------------------------------------------
 485    be % cv % size_jb = 0
 486    be % cv % size_je = 0
 487    be % cv % size_jp = 0
 488    be % cv % size_js = 0
 489    be % cv % size_jl = 0
 490    be % cv % size_jt = 0
 491
 492    !---------------------------------------------------------------------------
 493    ! [5.1] Set up background errors (be):
 494    !---------------------------------------------------------------------------
 495    if (use_background_errors .and. multi_inc /= 1) then
 496       call da_setup_background_errors (grid, be)
 497    else
 498       be%ncv_mz = 5
 499       allocate(be%cv_mz(5))
 500       be%cv_mz(:) = 0
 501       be % ne = ensdim_alpha
 502       be % v1 % mz = 0
 503       be % v2 % mz = 0
 504       be % v3 % mz = 0
 505       be % v4 % mz = 0
 506       be % v5 % mz = 0
 507       be % v6  % mz = 0
 508       be % v7  % mz = 0
 509       be % v8  % mz = 0
 510       be % v9  % mz = 0
 511       be % v10 % mz = 0
 512       be % v11 % mz = 0
 513    end if
 514
 515     ! overwrite variables defined in da_setup_cv.inc set in the call to da_setup_background_errors
 516     if ( anal_type_hybrid_dual_res ) then
 517        be % cv % size_alphac = (ite_int - its_int + 1) * (jte_int - jts_int + 1) * be % alpha % mz * be % ne
 518        be % cv % size_je = be % cv % size_alphac
 519        cv_size_domain_je = (ide_int - ids_int + 1) * (jde_int - jds_int + 1) * be % alpha % mz * be % ne
 520     endif
 521
 522    !---------------------------------------------------------------------------
 523    ! [5.2] Set up observation bias correction (VarBC):
 524    !---------------------------------------------------------------------------
 525 #if defined(RTTOV) || defined(CRTM)
 526    if (use_rad .and. (use_varbc.or.freeze_varbc)) call da_varbc_init(iv, be)
 527 #endif
 528
 529    if ( be%cv%size_jb == 0 .and. be%cv%size_je == 0 .and. be%cv%size_jp > 0 ) then
 530       offline_varbc = .true.
 531    else
 532       offline_varbc = .false.
 533    end if
 534
 535    if (use_tamdarobs .and. use_varbc_tamdar) then
 536         call da_varbc_tamdar_init(iv)
 537    else
 538         use_varbc_tamdar = .false.
 539    end if
 540    if (use_varbc_tamdar) call da_varbc_tamdar_pred(iv, be, ob)
 541
 542    !---------------------------------------------------------------------------
 543    ! [5.3] Set up satellite control variable:
 544    !---------------------------------------------------------------------------
 545 #if defined(RTTOV) || defined(CRTM)
 546    if (ANY(use_satcv)) call da_setup_satcv(iv, be)
 547 #endif
 548
 549    !---------------------------------------------------------------------------
 550    ! [5.4] Total control variable:
 551    !---------------------------------------------------------------------------
 552    be % cv % size = be%cv%size_jb + be%cv%size_je + be%cv%size_jp + be%cv%size_js + be%cv%size_jl + be%cv%size_jt
 553    cv_size = be % cv % size
 554
 555    !---------------------------------------------------------------------------
 556    ! [6.0] Set up ensemble perturbation input:
 557    !---------------------------------------------------------------------------
 558
 559       grid % ep % ne = be % ne
 560       if (use_background_errors .and. be % ne > 0) then
 561
 562          call date_and_time(date=this_date, time=this_time)
 563          write(unit=message(1),fmt='(a,a8,1x,a12)') &
 564             '   Begin reading ensemble perturbations at ', this_date, &
 565             this_time(1:2)//':'//this_time(3:4)//' '//this_time(5:10)
 566          call da_message(message(1:1))
 567
 568          select case ( ep_format )
 569             case ( 1 ) ! original behavior
 570                if ( ep_para_read == 0 ) then
 571                   ! each ep file is for one variable and one member
 572                   call da_setup_flow_predictors ( &
 573                      ide_ens, jde_ens, kde_ens, be % ne, grid%ep, &
 574                      its_int, ite_int, jts_int, jte_int, kts_int, kte_int )
 575                else if ( ep_para_read == 1 ) then
 576                   ! para_read_opt1, breaks total number of ep into reading bins
 577                   call da_setup_flow_predictors_para_read_opt1 ( &
 578                      ide_ens, jde_ens, kde_ens, be % ne, grid%ep, &
 579                      its_int, ite_int, jts_int, jte_int, kts_int, kte_int )
 580                end if
 581             case ( 11 )
 582                ! similar to ep_format=1 except data are in single precision
 583                ! each ep file is for one variable and one member
 584                call da_setup_flow_predictors ( &
 585                   ide_ens, jde_ens, kde_ens, be % ne, grid%ep, &
 586                   its_int, ite_int, jts_int, jte_int, kts_int, kte_int )
 587             case ( 2 )
 588                ! full-domain
 589                ! each ep file in single precision is for one variable and all members
 590                call da_setup_flow_predictors_ep_format2 ( &
 591                   ide_ens, jde_ens, kde_ens, be % ne, grid%ep, &
 592                   its_int, ite_int, jts_int, jte_int, kts_int, kte_int )
 593             case ( 3 )
 594                ! decomposed sub-domain
 595                ! each ep file in single precision is for one variable and all members
 596                call da_setup_flow_predictors_ep_format3 ( &
 597                   ide_ens, jde_ens, kde_ens, be % ne, grid%ep, &
 598                   its_int, ite_int, jts_int, jte_int, kts_int, kte_int )
 599          end select
 600
 601          call date_and_time(date=this_date, time=this_time)
 602          write(unit=message(1),fmt='(a,a8,1x,a12)') &
 603             '  Finish reading ensemble perturbations at ', this_date, &
 604             this_time(1:2)//':'//this_time(3:4)//' '//this_time(5:10)
 605          call da_message(message(1:1))
 606
 607       end if
 608
 609    !---------------------------------------------------------------------------
 610    ! [7.0] Setup control variable (cv):
 611    !---------------------------------------------------------------------------
 612
 613 !  Dynamically allocate the variables which don't rely on ntmax
 614       allocate (cvt(1:cv_size))
 615       allocate (xhat(1:cv_size))
 616 !      if (use_lanczos) then
 617 !        allocate (full_eignvec(cv_size))
 618 !      end if
 619
 620       !------------------------------------------------------
 621       ! set CV to random noise ("RANDOMCV")
 622       !------------------------------------------------------
 623       if (anal_type_randomcv) then
 624          it = 1
 625          ! Initialize random number generator and scalars
 626          call da_random_seed
 627          do i= 1, n_randomcv
 628             write(ci,'(i3.3)') i
 629             write(unit=message(1),fmt='(a,a)') &
 630                '  Setting randomcv for wrfvar_output_randomcv.e', trim(ci)
 631             call da_message(message(1:1))
 632             call da_set_randomcv    (cv_size, cvt)
 633             call da_transform_vtox  (grid,cv_size,xbx,be,grid%ep,cvt,grid%vv,grid%vp)
 634             call da_transform_xtoxa (grid)
 635             call da_transfer_xatoanalysis (it, xbx, grid, config_flags)
 636             call da_update_firstguess(grid,'wrfvar_output_randomcv.e'//trim(ci))
 637             if ( i < n_randomcv ) then
 638                ! restore the original grid info for the next realization of randomcv
 639                call da_med_initialdata_input( grid , config_flags, 'fg')
 640                call da_setup_firstguess(xbx, grid, config_flags, .false.)
 641             end if
 642          end do
 643          ! Done with randomcv.
 644          ! Set the following to skip some code to get to the deallocation part.
 645          max_ext_its = 0
 646       end if  !anal_type_randomcv
 647
 648 ! mri-4dvar: if multi_inc /= 2: run normal 3D/4D-Var
 649 !------------------------------------------------------------------------
 650     ! cvt is outer loop control variable, it is zero for the first outer loop,
 651     ! but non-zero from the second outer loop in normal 3d/4dvar.
 652     ! for MRI-4DVar, vp from the previous outer loop needs to be read in,
 653     ! then perform the inverse transform to derive cvt
 654     !-----------------------------------------------------
 655       call da_initialize_cv (cv_size, cvt)
 656       call da_zero_vp_type (grid%vp)
 657       call da_zero_vp_type (grid%vv)
 658
 659       if ( multi_inc == 2 ) then
 660          if ( max_ext_its > 1 ) then
 661            max_ext_its=1
 662            write(unit=message(1),fmt='(a)') "Re-set max_ext_its = 1 for multi_inc==2"
 663            call da_message(message(1:1))
 664          end if
 665
 666       ! read vp files for different PEs
 667       !----------------------------------
 668          write(unit=vpfile,fmt='(a,i4.4)') 'vp_input.',myproc
 669          inquire(file=trim(vpfile), exist=ex)
 670          if ( ex ) then
 671            call da_get_unit(vp_unit)
 672            open(unit=vp_unit,file=trim(vpfile),iostat=iost,form='UNFORMATTED',status='OLD')
 673            if (iost /= 0) then
 674              write(unit=message(1),fmt='(A,I5,A)') &
 675                 "Error ",iost," opening vp file "//trim(vpfile)
 676              call da_error(__FILE__,__LINE__,message(1:1))
 677            end if
 678            if ( use_interpolate_cvt ) then         ! works for CV3?, 3D RF
 679                write(unit=message(1),fmt='(a)') 'Reading vv from : '//trim(vpfile)
 680            elseif ( use_inverse_squarerootb ) then ! works for CV5,6,7, vertical EOF
 681                write(unit=message(1),fmt='(a)') 'Reading vp from : '//trim(vpfile)
 682            end if
 683            call da_message(message(1:1))
 684            read(vp_unit) i1, i2, i3, i4, i5, i6 ! read dimension of patch for current processor
 685            allocate( v1(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1) )
 686            allocate( v2(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1) )
 687            allocate( v3(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1) )
 688            allocate( v4(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1) )
 689            allocate( v5(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1) )
 690            if ( cloud_cv_options >= 2 ) then
 691              allocate( v6(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1) )
 692              allocate( v7(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1) )
 693              allocate( v8(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1) )
 694              allocate( v9(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1) )
 695              allocate( v10(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1) )
 696            end if
 697            if ( use_cv_w ) allocate( v11(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1) )
 698
 699            read(vp_unit) v1, v2, v3, v4, v5
 700            if ( cloud_cv_options >= 2 ) read(vp_unit) v6, v7, v8, v9, v10
 701            if ( use_cv_w ) read(vp_unit) v11
 702
 703            if ( use_interpolate_cvt ) then
 704              grid%vv%v1(ips:ipe,jps:jpe,kps:kpe) = v1(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
 705              grid%vv%v2(ips:ipe,jps:jpe,kps:kpe) = v2(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
 706              grid%vv%v3(ips:ipe,jps:jpe,kps:kpe) = v3(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
 707              grid%vv%v4(ips:ipe,jps:jpe,kps:kpe) = v4(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
 708              grid%vv%v5(ips:ipe,jps:jpe,kps:kpe) = v5(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
 709              if ( cloud_cv_options >= 2 ) then
 710                grid%vv%v6(ips:ipe,jps:jpe,kps:kpe) = v6(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
 711                grid%vv%v7(ips:ipe,jps:jpe,kps:kpe) = v7(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
 712                grid%vv%v8(ips:ipe,jps:jpe,kps:kpe) = v8(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
 713                grid%vv%v9(ips:ipe,jps:jpe,kps:kpe) = v9(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
 714                grid%vv%v10(ips:ipe,jps:jpe,kps:kpe) = v10(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
 715              end if
 716              if ( use_cv_w ) then
 717                grid%vv%v11(ips:ipe,jps:jpe,kps:kpe) = v11(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
 718              end if
 719              call da_vv_to_cv( grid%vv, grid%xp, be%cv_mz, be%ncv_mz, cv_size, cvt )
 720            elseif ( use_inverse_squarerootb ) then
 721              grid%vp%v1(ips:ipe,jps:jpe,kps:kpe) = v1(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
 722              grid%vp%v2(ips:ipe,jps:jpe,kps:kpe) = v2(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
 723              grid%vp%v3(ips:ipe,jps:jpe,kps:kpe) = v3(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
 724              grid%vp%v4(ips:ipe,jps:jpe,kps:kpe) = v4(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
 725              grid%vp%v5(ips:ipe,jps:jpe,kps:kpe) = v5(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
 726              if ( cloud_cv_options >= 2 ) then
 727                grid%vp%v6(ips:ipe,jps:jpe,kps:kpe) = v6(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
 728                grid%vp%v7(ips:ipe,jps:jpe,kps:kpe) = v7(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
 729                grid%vp%v8(ips:ipe,jps:jpe,kps:kpe) = v8(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
 730                grid%vp%v9(ips:ipe,jps:jpe,kps:kpe) = v9(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
 731                grid%vp%v10(ips:ipe,jps:jpe,kps:kpe) = v10(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
 732              end if
 733              if ( use_cv_w ) then ! vertical stagging +1?
 734                grid%vp%v11(ips:ipe,jps:jpe,kps:kpe) = v11(1:i2-i1+1, 1:i4-i3+1, 1:i6-i5+1)
 735              end if
 736              !call da_write_vp(grid,grid%vp,'vp_input.global ') ! to verify correctness
 737              print '(/10X,"===> Use inverse transform of square-root B for outer-loop=",i2)', it
 738              if ( cv_options == 3 ) then
 739                 write(unit=message(1),fmt='(A,I5,A)') &
 740                    "Error: inverse U transform not for cv_options = 3"
 741                 call da_error(__FILE__,__LINE__,message(1:1))
 742              end if
 743              call da_transform_vtox_inv (grid,be%cv%size_jb,xbx,be,grid%ep,cvt(1:be%cv%size_jb),grid%vv,grid%vp)
 744            end if
 745
 746            deallocate( v1 )
 747            deallocate( v2 )
 748            deallocate( v3 )
 749            deallocate( v4 )
 750            deallocate( v5 )
 751            if ( cloud_cv_options >= 2 ) then
 752              deallocate( v6 )
 753              deallocate( v7 )
 754              deallocate( v8 )
 755              deallocate( v9 )
 756              deallocate( v10 )
 757            end if
 758            if ( use_cv_w ) deallocate( v11 )
 759
 760            close(vp_unit)
 761            call da_free_unit(vp_unit)
 762
 763          else
 764            write(unit=message(1),fmt='(a)') "vp files '"//trim(vpfile)//"' does not exists, initiallizing cvt."
 765            call da_message(message(1:1))
 766            call da_initialize_cv (cv_size, cvt) ! perhaps better use da_error
 767          end if
 768
 769          call da_zero_vp_type (grid%vv)
 770          call da_zero_vp_type (grid%vp)
 771
 772       end if
 773 ! mri-4dvar -------------------------------------------
 774
 775       if ( var4d ) then
 776 #ifdef VAR4D
 777          call da_zero_vp_type (grid%vv6)
 778          call da_zero_vp_type (grid%vp6)
 779 #endif
 780       end if
 781
 782    !---------------------------------------------------------------------------
 783    ! [8] Outerloop
 784    !---------------------------------------------------------------------------
 785
 786    j_grad_norm_target = 1.0
 787    do it = 1, max_ext_its
 788
 789 !  Dynamically allocate the variables which depend on ntmax
 790       if (use_lanczos) then
 791          allocate (qhat(1:cv_size, 0:ntmax(it)))
 792          allocate (eignvec(ntmax(it), ntmax(it)))
 793          allocate (eignval(ntmax(it)))
 794       end if
 795
 796 ! Re-scale the variances and the scale-length for outer-loop > 1:
 797       if (it > 1 .and. (cv_options == 5 .or. cv_options == 7)) then
 798          print '(/10X,"===> Re-set BE SCALINGS for outer-loop=",i2)', it
 799          call  da_scale_background_errors ( be, it )
 800       else if (it > 1 .and. cv_options == 3) then
 801          print '(/10X,"===> Re-set CV3 BE SCALINGS for outer-loop=",i2)', it
 802          call  da_scale_background_errors_cv3 ( grid, be, it )
 803       endif
 804
 805       call da_initialize_cv (cv_size, xhat)
 806
 807 ! mri-4dvar----------------------
 808 ! Apply inverse transform of squareroot(B) for full-resolution non-stop Var
 809 ! from 2nd outer loop, this is to check correctness of inverse U transform
 810 ! does not apply this setting for real world application
 811 !-----------------------------
 812       if (multi_inc == 0 .and. it > 1 .and. use_inverse_squarerootb .and. cv_options /= 3) then
 813          print '(/10X,"===> Use inverse transform of square-root B for outer-loop=",i2)', it
 814          call da_transform_vtox_inv (grid,be%cv%size_jb,xbx,be,grid%ep,cvt(1:be%cv%size_jb),grid%vv,grid%vp)
 815       endif
 816
 817 ! Reinitialize cvt=0 for full-resolution non-stop Var for each loop
 818 !------another option not tested --------------
 819       if (multi_inc == 0 .and. it > 1 .and. use_interpolate_cvt) then
 820          print '(/10X,"===> Reinitialize cvt as zeros for outer loop ",i2)', it
 821          call da_initialize_cv (cv_size, cvt)
 822       endif
 823 ! mri-4dvar------------------------
 824
 825       ! [8.1] Calculate nonlinear model trajectory
 826
 827 !     if (var4d .and. multi_inc /= 2 ) then
 828       if (var4d) then
 829 #ifdef VAR4D
 830          if (it > 1) then
 831             call kj_swap (grid%u_2, model_grid%u_2, &
 832                           grid%xp%ims, grid%xp%ime, grid%xp%jms, grid%xp%jme, grid%xp%kms, grid%xp%kme)
 833             call kj_swap (grid%v_2, model_grid%v_2, &
 834                           grid%xp%ims, grid%xp%ime, grid%xp%jms, grid%xp%jme, grid%xp%kms, grid%xp%kme)
 835             call kj_swap (grid%w_2, model_grid%w_2, &
 836                           grid%xp%ims, grid%xp%ime, grid%xp%jms, grid%xp%jme, grid%xp%kms, grid%xp%kme)
 837             call kj_swap (grid%t_2, model_grid%t_2, &
 838                           grid%xp%ims, grid%xp%ime, grid%xp%jms, grid%xp%jme, grid%xp%kms, grid%xp%kme)
 839             call kj_swap (grid%ph_2, model_grid%ph_2, &
 840                           grid%xp%ims, grid%xp%ime, grid%xp%jms, grid%xp%jme, grid%xp%kms, grid%xp%kme)
 841             call kj_swap (grid%p, model_grid%p, &
 842                           grid%xp%ims, grid%xp%ime, grid%xp%jms, grid%xp%jme, grid%xp%kms, grid%xp%kme)
 843             model_grid%mu_2 = grid%mu_2
 844             model_grid%t2 = grid%t2
 845             model_grid%th2 = grid%th2
 846             model_grid%q2 = grid%q2
 847             model_grid%u10 = grid%u10
 848             model_grid%v10 = grid%v10
 849             model_grid%tsk = grid%tsk
 850             model_grid%psfc = grid%psfc
 851
 852             do i = PARAM_FIRST_SCALAR, num_moist
 853                call kj_swap (grid%moist(:,:,:,i), model_grid%moist(:,:,:,i), &
 854                              grid%xp%ims, grid%xp%ime, grid%xp%jms, grid%xp%jme, grid%xp%kms, grid%xp%kme)
 855             enddo
 856
 857             ! Update boundary condition of model
 858
 859             model_grid%u_bxs = grid%u_bxs
 860             model_grid%u_bxe = grid%u_bxe
 861             model_grid%u_bys = grid%u_bys
 862             model_grid%u_bye = grid%u_bye
 863             model_grid%v_bxs = grid%v_bxs
 864             model_grid%v_bxe = grid%v_bxe
 865             model_grid%v_bys = grid%v_bys
 866             model_grid%v_bye = grid%v_bye
 867             model_grid%w_bxs = grid%w_bxs
 868             model_grid%w_bxe = grid%w_bxe
 869             model_grid%w_bys = grid%w_bys
 870             model_grid%w_bye = grid%w_bye
 871             model_grid%t_bxs = grid%t_bxs
 872             model_grid%t_bxe = grid%t_bxe
 873             model_grid%t_bys = grid%t_bys
 874             model_grid%t_bye = grid%t_bye
 875             model_grid%mu_bxs = grid%mu_bxs
 876             model_grid%mu_bxe = grid%mu_bxe
 877             model_grid%mu_bys = grid%mu_bys
 878             model_grid%mu_bye = grid%mu_bye
 879             model_grid%ph_bxs = grid%ph_bxs
 880             model_grid%ph_bxe = grid%ph_bxe
 881             model_grid%ph_bys = grid%ph_bys
 882             model_grid%ph_bye = grid%ph_bye
 883             model_grid%moist_bxs = grid%moist_bxs
 884             model_grid%moist_bxe = grid%moist_bxe
 885             model_grid%moist_bys = grid%moist_bys
 886             model_grid%moist_bye = grid%moist_bye
 887
 888             model_grid%u_btxs = grid%u_btxs
 889             model_grid%u_btxe = grid%u_btxe
 890             model_grid%u_btys = grid%u_btys
 891             model_grid%u_btye = grid%u_btye
 892             model_grid%v_btxs = grid%v_btxs
 893             model_grid%v_btxe = grid%v_btxe
 894             model_grid%v_btys = grid%v_btys
 895             model_grid%v_btye = grid%v_btye
 896             model_grid%w_btxs = grid%w_btxs
 897             model_grid%w_btxe = grid%w_btxe
 898             model_grid%w_btys = grid%w_btys
 899             model_grid%w_btye = grid%w_btye
 900             model_grid%t_btxs = grid%t_btxs
 901             model_grid%t_btxe = grid%t_btxe
 902             model_grid%t_btys = grid%t_btys
 903             model_grid%t_btye = grid%t_btye
 904             model_grid%mu_btxs = grid%mu_btxs
 905             model_grid%mu_btxe = grid%mu_btxe
 906             model_grid%mu_btys = grid%mu_btys
 907             model_grid%mu_btye = grid%mu_btye
 908             model_grid%ph_btxs = grid%ph_btxs
 909             model_grid%ph_btxe = grid%ph_btxe
 910             model_grid%ph_btys = grid%ph_btys
 911             model_grid%ph_btye = grid%ph_btye
 912             model_grid%moist_btxs = grid%moist_btxs
 913             model_grid%moist_btxe = grid%moist_btxe
 914             model_grid%moist_btys = grid%moist_btys
 915             model_grid%moist_btye = grid%moist_btye
 916
 917             ! Turn off model boundary reading as we already provide a new one.
 918             call da_model_lbc_off
 919          endif
 920
 921          call nl_set_var4d_run (head_grid%id, .true.)
 922          call da_nl_model(it)
 923 !      elseif (var4d .and. multi_inc == 2 ) then
 924 #else
 925          write(unit=message(1),fmt='(A)')'Please re-compile the code with 4dvar option'
 926          call da_error(__FILE__,__LINE__,message(1:1))
 927 #endif
 928       end if
 929
 930       ! [8.2] Calculate innovation vector (O-B):
 931
 932       num_qcstat_conv=0
 933       call da_get_innov_vector (it, num_qcstat_conv, ob, iv, grid , config_flags)
 934       if ( multi_inc == 1 ) then
 935          if (trace_use) call da_trace_exit ("da_solve")
 936          return
 937       end if
 938
 939       if (test_transforms .or. test_gradient) then
 940
 941          if (test_gradient) then
 942             call da_allocate_y (iv, re)
 943             call da_allocate_y (iv, y)
 944             call da_check_gradient (grid, config_flags, cv_size, xhat, cvt, 1.0e-10, 8, &
 945                  xbx, be, iv, y, re, j)
 946             call da_deallocate_y (re)
 947             call da_deallocate_y (y)
 948          endif
 949
 950          if (test_transforms) then
 951             call da_check (grid, config_flags, cv_size, xbx, be, grid%ep, iv, &
 952                         grid%vv, grid%vp, y)
 953          endif
 954
 955          if (trace_use) call da_trace_exit("da_solve")
 956          return
 957
 958       end if
 959
 960       ! [8.4] Minimize cost function:
 961
 962       call da_allocate_y (iv, re)
 963       call da_allocate_y (iv, y)
 964
 965       if (use_lanczos) then
 966          if (read_lanczos) then
 967             call da_lanczos_io('r',cv_size,ntmax(it),neign,eignvec,eignval,qhat)
 968
 969             call da_kmat_mul(grid,config_flags,it,cv_size,xbx, &
 970                              be,iv,xhat,qhat,cvt,re,y,j,eignvec,eignval,neign)
 971           ! Output Cost Function
 972             call da_calculate_j(it, 1, cv_size, be%cv%size_jb, be%cv%size_je, be%cv%size_jp, &
 973                              be%cv%size_jl, be%cv%size_jt, xbx, be, iv, xhat, cvt, re, y, j, grid, config_flags )
 974
 975          else
 976             call da_minimise_lz(grid, config_flags, it, cv_size, xbx,&
 977                              be, iv, j_grad_norm_target, xhat, qhat, cvt, re, y, j, eignvec, eignval, neign )
 978
 979
 980          end if
 981
 982          if (write_lanczos) call da_lanczos_io('w',cv_size,ntmax(it),neign,eignvec,eignval,qhat)
 983
 984          if (adj_sens) call da_sensitivity(grid,config_flags,it,cv_size,xbx, &
 985                                            be,iv,xhat,qhat,cvt,y,eignvec,eignval,neign )
 986
 987
 988       else
 989
 990          call da_minimise_cg( grid, config_flags, it, be % cv % size, &
 991               xbx, be, iv, j_grad_norm_target, xhat, cvt, re, y, j)
 992       end if
 993
 994       ! Update outer-loop control variable
 995       cvt = cvt + xhat
 996       if ( multi_inc == 2 .and. use_interpolate_cvt ) then ! obsolete option
 997         call da_cv_to_vv( cv_size, cvt, be%cv_mz, be%ncv_mz, grid%vv )
 998         call da_write_vp(grid,grid%vv,'vp_output.global') ! wrtie vv to vp file
 999       end if
1000       !------------------------------------------------------------------------
1001
1002       ! [8.5] Update latest analysis solution:
1003
1004       if (.not. var4d) then
1005          ! as of V3.9, vp%alpha is not included in the vptox (called by vtox) calculation,
1006          ! here grid%xa contains only the static increment
1007
1008 #if (WRF_CHEM == 1)
1009    if (cv_options == 7 .and. chem_cv_options>=10) then
1010          call da_transform_vtox (grid,cv_size,xbx,be,grid%ep,xhat,grid%vv,grid%vp,vchem=grid%vchem)
1011    end if
1012 #else
1013          call da_transform_vtox (grid,cv_size,xbx,be,grid%ep,xhat,grid%vv,grid%vp)
1014 #endif
1015
1016          ! adding the ensemble contribution
1017          if (be % ne > 0 .and. alphacv_method == alphacv_method_xa) then
1018             if ( use_4denvar ) then
1019                ! add it to the analysis time ifgat_ana
1020                ! ifgat_ana is declared in da_control.f90 and defined in da_get_time_slots.inc
1021                write(unit=message(1),fmt='(a)')    ''
1022                write(unit=message(2),fmt='(a,i3)') 'Applying ensemble contribution to FGAT time: ', ifgat_ana
1023                call da_message(message(1:2))
1024                call da_transform_vpatox (grid,be%ne,grid%ep,grid%vp,ifgat_ana)
1025             else
1026                ! ep does not have time dimension for non-4denvar
1027                call da_transform_vpatox (grid,be%ne,grid%ep,grid%vp)
1028             end if
1029             call da_add_xa(grid%xa, grid%xa_ens) !grid%xa = grid%xa + xa_ens
1030          end if !ne>0
1031       else
1032          call da_transform_vtox (grid,be%cv%size_jb,xbx,be,grid%ep,xhat(1:be%cv%size_jb),grid%vv,grid%vp)
1033          call da_transform_vpatox (grid,be%ne,grid%ep,grid%vp)
1034       endif
1035
1036 ! mri-4dvar--------------------------
1037       if (multi_inc == 2 .and. use_inverse_squarerootb) then
1038          call da_write_vp(grid,grid%vp,'vp_output.global') ! write vp to vp file
1039       end if
1040 ! mri-4dvar--------------------------
1041
1042       call da_transform_xtoxa (grid)
1043
1044       ! [8.6] Only when use_radarobs = .false. and calc_w_increment =.true.,
1045       !       the w_increment need to be diagnosed:
1046
1047       if (calc_w_increment .and. .not. use_radarobs .and. .not. var4d) then
1048          call da_uvprho_to_w_lin (grid)
1049
1050 #ifdef DM_PARALLEL
1051 #include "HALO_RADAR_XA_W.inc"
1052 #endif
1053       end if
1054
1055       ! [8.7] Write out diagnostics
1056
1057       call da_write_diagnostics (it, grid, num_qcstat_conv, ob, iv, re, y, j)
1058
1059       ! Write "clean" QCed observations if requested:
1060       if (anal_type_qcobs) then
1061          ! if (it == 1) then
1062           if (write_mod_filtered_obs) then
1063             call da_write_modified_filtered_obs (grid, ob, iv, &
1064                coarse_ix, coarse_jy, start_x, start_y)
1065           else
1066             call da_write_filtered_obs (it, grid, ob, iv, &
1067                coarse_ix, coarse_jy, start_x, start_y)
1068           end if
1069          ! end if
1070       end if
1071
1072       ! [8.7.1] Write Ascii radar OMB and OMA file
1073
1074       if ( use_radarobs .and. write_oa_radar_ascii ) then
1075          call da_write_oa_radar_ascii (ob,iv,re,it)
1076       end if
1077
1078       ! [8.3] Interpolate x_g to low resolution grid
1079
1080       ! [8.8] Write Ascii radiance OMB and OMA file
1081
1082 #if defined(CRTM) || defined(RTTOV)
1083       if (use_rad .and. write_oa_rad_ascii) then
1084          call da_write_oa_rad_ascii (it,ob,iv,re)
1085       end if
1086 #endif
1087
1088       ! [8.9] Update VarBC parameters and write output file
1089 #if defined(CRTM) || defined(RTTOV)
1090       if ( use_rad .and. (use_varbc.or.freeze_varbc) ) &
1091                 call da_varbc_update(it, cv_size, xhat, iv)
1092 #endif
1093
1094       if (use_varbc_tamdar) &
1095          call da_varbc_tamdar_update(cv_size, xhat, iv)
1096
1097       !------------------------------------------------------------------------
1098       ! [8.10] Output WRFVAR analysis and analysis increments:
1099       !------------------------------------------------------------------------
1100
1101       if ( .not. offline_varbc ) &
1102       call da_transfer_xatoanalysis (it, xbx, grid, config_flags)
1103
1104      if ( it < max_ext_its .and. print_detail_outerloop ) then
1105         write(outerloop,'(i2.2)') it
1106         call da_update_firstguess(grid,'wrfvar_output_'//outerloop)
1107 #ifdef VAR4D
1108         !if (var4d) call da_med_initialdata_output_lbc (grid , config_flags, 'wrfvar_bdyout_'//outerloop)
1109 #endif
1110      end if
1111
1112      call da_deallocate_y (re)
1113      call da_deallocate_y (y)
1114
1115
1116    ! Deallocate arrays which depend on ntmax
1117      if (use_lanczos) then
1118         deallocate (qhat)
1119         deallocate (eignvec)
1120         deallocate (eignval)
1121      end if
1122
1123
1124    end do
1125
1126    ! output wrfvar analysis
1127
1128    if ((config_flags%real_data_init_type == 1) .or. &
1129        (config_flags%real_data_init_type == 3)) then
1130       if ( .not. offline_varbc ) &
1131       call da_update_firstguess(input_grid)
1132 #ifdef VAR4D
1133       !if (var4d) call da_med_initialdata_output_lbc (head_grid , config_flags)
1134       if ( var4d_lbc ) then
1135          call domain_clock_get (grid, stop_timestr=timestr1)
1136          call domain_clock_set( grid, current_timestr=timestr1 )
1137          call da_med_initialdata_input (grid, config_flags, 'fg02')
1138          call da_setup_firstguess(xbx, grid, config_flags, .false. )
1139          shuffle = grid%xa
1140          jl_start    = be%cv%size_jb + be%cv%size_je + be%cv%size_jp + 1
1141          jl_end      = be%cv%size_jb + be%cv%size_je + be%cv%size_jp + be%cv%size_jl
1142          grid%xa  = grid%x6a
1143          call da_transform_vtox(grid, be%cv%size_jl, xbx, be, grid%ep, &
1144               xhat(jl_start:jl_end), grid%vv6, grid%vp6)
1145          grid%xa  = shuffle
1146          call da_transfer_xatoanalysis (it, xbx, grid, config_flags)
1147          call da_update_firstguess (grid, 'ana02')
1148          call domain_clock_get (grid, start_timestr=timestr1)
1149          call domain_clock_set( grid, current_timestr=timestr1 )
1150       endif
1151 #endif
1152       call med_shutdown_io (input_grid, config_flags)
1153    end if
1154
1155    !---------------------------------------------------------------------------
1156    ! [9.0] Tidy up:
1157    !---------------------------------------------------------------------------
1158
1159    deallocate (cvt)
1160    deallocate (xhat)
1161 !   if (use_lanczos) then
1162 !      deallocate (full_eignvec)
1163 !   end if
1164
1165    ! clean up radiance related arrays
1166 #if defined(RTTOV) || defined(CRTM)
1167    if (use_rad) then
1168       call da_deallocate_radiance (ob, iv, j)
1169       deallocate (time_slots)
1170 #ifdef RTTOV
1171       if (rtm_option == rtm_option_rttov) then
1172          deallocate (coefs)
1173          deallocate (opts)
1174       end if
1175 #endif
1176    end if
1177 #endif
1178
1179    if (var4d .and. use_rainobs) deallocate(fgat_rain_flags)
1180    call da_deallocate_observations (iv)
1181    call da_deallocate_y (ob)
1182    if (use_background_errors .and. jb_factor>0.0) call da_deallocate_background_errors (be)
1183
1184    if (xbx%pad_num > 0) then
1185       deallocate (xbx%pad_loc)
1186       deallocate (xbx%pad_pos)
1187    end if
1188
1189    deallocate (xbx % fft_factors_x)
1190    deallocate (xbx % fft_factors_y)
1191    deallocate (xbx % fft_coeffs)
1192    deallocate (xbx % trig_functs_x)
1193    deallocate (xbx % trig_functs_y)
1194
1195    if (global) then
1196       deallocate (xbx%coslat)
1197       deallocate (xbx%sinlat)
1198       deallocate (xbx%coslon)
1199       deallocate (xbx%sinlon)
1200       deallocate (xbx%int_wgts)
1201       deallocate (xbx%alp)
1202       deallocate (xbx%wsave)
1203       if (jts == jds) then
1204          deallocate (cos_xls)
1205          deallocate (sin_xls)
1206       end if
1207
1208       if (jte == jde) then
1209          deallocate (cos_xle)
1210          deallocate (sin_xle)
1211       end if
1212    end if
1213
1214    if ( anal_type_hybrid_dual_res ) deallocate(aens_locs)
1215
1216    if (use_varbc_tamdar) then
1217       deallocate (iv%varbc_tamdar%nobs)
1218       deallocate (iv%varbc_tamdar%nobs_sum)
1219       deallocate (iv%varbc_tamdar%tail_id)
1220       deallocate (iv%varbc_tamdar%obs_sn)
1221       deallocate (iv%varbc_tamdar%ifuse)
1222       deallocate (iv%varbc_tamdar%index)
1223       deallocate (iv%varbc_tamdar%pred)
1224       deallocate (iv%varbc_tamdar%param)
1225       deallocate (iv%varbc_tamdar%bgerr)
1226       deallocate (iv%varbc_tamdar%vtox)
1227    end if
1228
1229    deallocate ( c1f )
1230    deallocate ( c2f )
1231    deallocate ( c3f )
1232    deallocate ( c4f )
1233    deallocate ( c1h )
1234    deallocate ( c2h )
1235    deallocate ( c3h )
1236    deallocate ( c4h )
1237
1238 #ifdef DM_PARALLEL
1239    call mpi_barrier (comm,ierr)
1240 #endif
1241
1242    if (trace_use) call da_trace_exit ("da_solve")
1243
1244 end subroutine da_solve
1245