a collection of WRFDA fixes related to ensdim_alpha use cases

var/da/da_main/da_solve.inc

@@ -59,6 +59,7 @@
    character(len=13) :: vpfile ! vp_input.0001
    integer   :: i1,i2,i3,i4,i5,i6
    logical :: ex
+   logical :: offline_varbc
 
    character(len=10) :: this_time
    character(len=8)  :: this_date
@@ -494,6 +495,9 @@
    if (use_background_errors .and. multi_inc /= 1) then
       call da_setup_background_errors (grid, be)
    else
+      be%ncv_mz = 5
+      allocate(be%cv_mz(5))
+      be%cv_mz(:) = 0
       be % ne = ensdim_alpha
       be % v1 % mz = 0
       be % v2 % mz = 0
@@ -522,6 +526,12 @@
    if (use_rad .and. (use_varbc.or.freeze_varbc)) call da_varbc_init(iv, be)
 #endif
 
+   if ( be%cv%size_jb == 0 .and. be%cv%size_je == 0 .and. be%cv%size_jp > 0 ) then
+      offline_varbc = .true.
+   else
+      offline_varbc = .false.
+   end if
+
    if (use_tamdarobs .and. use_varbc_tamdar) then
         call da_varbc_tamdar_init(iv)
    else
@@ -590,7 +600,7 @@
 
          call date_and_time(date=this_date, time=this_time)
          write(unit=message(1),fmt='(a,a8,1x,a12)') &
-            '  Fihish reading ensemble perturbations at ', this_date, &
+            '  Finish reading ensemble perturbations at ', this_date, &
             this_time(1:2)//':'//this_time(3:4)//' '//this_time(5:10)
          call da_message(message(1:1))
 
@@ -1005,17 +1015,22 @@
 
          ! adding the ensemble contribution
          if (be % ne > 0 .and. alphacv_method == alphacv_method_xa) then
-            ! add it to the analysis time ifgat_ana
-            ! ifgat_ana is declared in da_control.f90 and defined in da_get_time_slots.inc
-            write(unit=message(1),fmt='(a)')    ''
-            write(unit=message(2),fmt='(a,i3)') 'Applying ensemble contribution to FGAT time: ', ifgat_ana
-            call da_message(message(1:2))
-            call da_transform_vpatox (grid,be,grid%ep,grid%vp,ifgat_ana)
+            if ( use_4denvar ) then
+               ! add it to the analysis time ifgat_ana
+               ! ifgat_ana is declared in da_control.f90 and defined in da_get_time_slots.inc
+               write(unit=message(1),fmt='(a)')    ''
+               write(unit=message(2),fmt='(a,i3)') 'Applying ensemble contribution to FGAT time: ', ifgat_ana
+               call da_message(message(1:2))
+               call da_transform_vpatox (grid,be%ne,grid%ep,grid%vp,ifgat_ana)
+            else
+               ! ep does not have time dimension for non-4denvar
+               call da_transform_vpatox (grid,be%ne,grid%ep,grid%vp)
+            end if
             call da_add_xa(grid%xa, grid%xa_ens) !grid%xa = grid%xa + xa_ens
          end if !ne>0
       else
          call da_transform_vtox (grid,be%cv%size_jb,xbx,be,grid%ep,xhat(1:be%cv%size_jb),grid%vv,grid%vp)
-         call da_transform_vpatox (grid,be,grid%ep,grid%vp)
+         call da_transform_vpatox (grid,be%ne,grid%ep,grid%vp)
       endif
 
 ! mri-4dvar--------------------------
@@ -1083,6 +1098,7 @@
       ! [8.10] Output WRFVAR analysis and analysis increments:
       !------------------------------------------------------------------------
 
+      if ( .not. offline_varbc ) &
       call da_transfer_xatoanalysis (it, xbx, grid, config_flags)
 
      if ( it < max_ext_its .and. print_detail_outerloop ) then
@@ -1111,6 +1127,7 @@
 
    if ((config_flags%real_data_init_type == 1) .or. &
        (config_flags%real_data_init_type == 3)) then
+      if ( .not. offline_varbc ) &
       call da_update_firstguess(input_grid)
 #ifdef VAR4D
       !if (var4d) call da_med_initialdata_output_lbc (head_grid , config_flags)
@@ -1162,7 +1179,7 @@
    if (var4d .and. use_rainobs) deallocate(fgat_rain_flags)
    call da_deallocate_observations (iv)
    call da_deallocate_y (ob)
-   if (use_background_errors) call da_deallocate_background_errors (be)
+   if (use_background_errors .and. jb_factor>0.0) call da_deallocate_background_errors (be)
 
    if (xbx%pad_num > 0) then
       deallocate (xbx%pad_loc)

var/da/da_minimisation/da_calculate_gradj.inc

@@ -222,7 +222,7 @@ subroutine da_calculate_gradj(it, iter, cv_size, cv_size_jb, cv_size_je, cv_size
    if (use_divc) then
       call da_transform_vtox(grid, cv_size, xbx, be, grid%ep, cv, grid%vv, grid%vp)
       if ( be%ne > 0 .and. alphacv_method == alphacv_method_xa ) then
-         call da_transform_vpatox(grid, be, grid%ep, grid%vp)
+         call da_transform_vpatox(grid, be%ne, grid%ep, grid%vp)
          call da_add_xa(grid%xa, grid%xa_ens) !grid%xa = grid%xa + xa_ens
       end if
       call da_transform_xtoxa(grid)
@@ -232,7 +232,7 @@ subroutine da_calculate_gradj(it, iter, cv_size, cv_size_jb, cv_size_je, cv_size
       call da_divergence_constraint_adj(grid, grid%xa%u, grid%xa%v, inc_div)
       call da_transform_xtoxa_adj(grid)
       if (be % ne > 0 .and. alphacv_method == alphacv_method_xa) then
-         call da_transform_vpatox_adj(grid, be, grid%ep, grid%vp)
+         call da_transform_vpatox_adj(grid, be%ne, grid%ep, grid%vp)
       end if
       call da_zero_vp_type(grid%vp)
       call da_transform_vtox_adj(grid, cv_size, xbx, be, grid%ep, grid%vp, grid%vv, grad_jm)

var/da/da_minimisation/da_calculate_j.inc

@@ -324,7 +324,7 @@ subroutine da_calculate_j(it, iter, cv_size, cv_size_jb, cv_size_je, cv_size_jp,
    if (use_divc) then
       call da_transform_vtox(grid, cv_size, xbx, be, grid%ep, xhat+cv, grid%vv, grid%vp)
       if ( be%ne > 0 .and. alphacv_method == alphacv_method_xa ) then
-         call da_transform_vpatox(grid, be, grid%ep, grid%vp)
+         call da_transform_vpatox(grid, be%ne, grid%ep, grid%vp)
          call da_add_xa(grid%xa, grid%xa_ens) !grid%xa = grid%xa + xa_ens
       end if
       call da_transform_xtoxa(grid)

var/da/da_minimisation/da_minimise_cg.inc

@@ -51,6 +51,7 @@ subroutine da_minimise_cg(grid, config_flags,            &
    ! Variables for Conjugate Gradient preconditioning
    real                              :: precon(1:cv_size)      ! cv copy.
    real                              :: g_total, g_partial, jo_partial                          
+   integer                           :: nv_precon
    integer                           :: i, ii, nv, nn, istart, iend
    integer, allocatable              ::  sz(:)
 
@@ -62,22 +63,6 @@ subroutine da_minimise_cg(grid, config_flags,            &
    !-------------------------------------------------------------------------
    ! [1.0] Initialization:
    !-------------------------------------------------------------------------
-   allocate (sz(be%ncv_mz-2))
-   sz(1) = be%cv%size1
-   sz(2) = be%cv%size2
-   sz(3) = be%cv%size3
-   sz(4) = be%cv%size4
-   sz(5) = be%cv%size5
-   if ( be%ncv_mz > 10 ) then
-      sz(6)  = be%cv%size6
-      sz(7)  = be%cv%size7
-      sz(8)  = be%cv%size8
-      sz(9)  = be%cv%size9
-      sz(10) = be%cv%size10
-   end if
-   if ( be%ncv_mz == 8 .or. be%ncv_mz == 13 ) then
-      sz(be%ncv_mz-2) = be%cv%size11i
-   end if
 
    jp_start   = be % cv % size_jb + be % cv % size_je + 1
    jp_end     = be % cv % size_jb + be % cv % size_je + be % cv % size_jp
@@ -102,20 +87,46 @@ subroutine da_minimise_cg(grid, config_flags,            &
    if (precondition_cg) then
       g_total = da_dot(cv_size,ghat,ghat)
 
+      !cv_mz: v1, v2, v3, v4, v5, [v6, v7, v8, v9, v10], [v11], [alpha_cv, ne]
+      if ( be%ne > 0 ) then
+         nv_precon = be%ncv_mz - 1
+      else
+         nv_precon = be%ncv_mz - 2
+      end if
+      allocate (sz(nv_precon))
+      sz(1) = be%cv%size1
+      sz(2) = be%cv%size2
+      sz(3) = be%cv%size3
+      sz(4) = be%cv%size4
+      sz(5) = be%cv%size5
+      if ( be%ncv_mz > 10 ) then
+         sz(6)  = be%cv%size6
+         sz(7)  = be%cv%size7
+         sz(8)  = be%cv%size8
+         sz(9)  = be%cv%size9
+         sz(10) = be%cv%size10
+      end if
+      if ( be%ncv_mz == 8 .or. be%ncv_mz == 13 ) then ! has cv_w
+         sz(be%ncv_mz-2) = be%cv%size11i
+      end if
+      if ( be%ne > 0 ) then
+         sz(nv_precon) = be%cv%size_alphac
+      end if
+
       iend    = 0
-      do nv = 1, (be%ncv_mz - 2)
+      do nv = 1, nv_precon
          nn = sz(nv) / be%cv_mz(nv)
          do ii = 1, be%cv_mz(nv)
             istart     = iend + 1
             iend       = istart + nn - 1
             g_partial  = da_dot(nn, ghat(istart:iend), ghat(istart:iend))
-            jo_partial = j0_total / SUM(be%cv_mz(1:(be%ncv_mz-2)))
-            precon(istart:iend)=  1 / &
-               (1 + precondition_factor*(g_partial/g_total)/(jo_partial/j0_total))
+            jo_partial = j0_total / SUM(be%cv_mz(1:nv_precon))
+            precon(istart:iend)=  1.0 / &
+               (1.0 + precondition_factor*(g_partial/g_total)/(jo_partial/j0_total))
          end do
       end do
+      deallocate(sz)
    end if
-   deallocate(sz)
 
    phat  = - precon * ghat
 
@@ -168,7 +179,7 @@ subroutine da_minimise_cg(grid, config_flags,            &
     ! Orthonormalize new gradient (using modified Gramm-Schmidt algorithm)
       if (orthonorm_gradient) then
          do i = iter-1, 0, -1
-            gdot = da_dot_cv(cv_size, ghat, qhat(i)%values, grid, be%cv_mz, be%ncv_mz &
+            gdot = da_dot_cv(cv_size, precon*ghat, qhat(i)%values, grid, be%cv_mz, be%ncv_mz &
 #if (WRF_CHEM == 1)
                ,be%cv_mz_chemic &
 #endif

var/da/da_minimisation/da_transform_vtod_wpec.inc

@@ -21,7 +21,7 @@ subroutine da_transform_vtod_wpec(cv_size, be, ep, cv, vp, vv, xbx, grid)
       call da_transform_vtox(grid, cv_size, xbx, be, ep, cv, vv, vp)
 
       if ( be%ne > 0 .and. alphacv_method == alphacv_method_xa ) then
-         call da_transform_vpatox(grid, be, ep, vp)
+         call da_transform_vpatox(grid, be%ne, ep, vp)
          call da_add_xa(grid%xa, grid%xa_ens) !grid%xa = grid%xa + xa_ens
       end if
 

var/da/da_minimisation/da_transform_vtod_wpec_adj.inc

@@ -27,7 +27,7 @@ subroutine da_transform_vtod_wpec_adj(cv_size, be, ep, cv, vp, vv, xbx, grid)
       call da_transform_xtoxa_adj(grid, .true.)
 
       if (be % ne > 0 .and. alphacv_method == alphacv_method_xa) then
-         call da_transform_vpatox_adj(grid, be, ep, vp)
+         call da_transform_vpatox_adj(grid, be%ne, ep, vp)
       end if
 
       call da_transform_vtox_adj(grid, cv_size, xbx, be, ep, vp, vv, cv)

var/da/da_minimisation/da_transform_vtoy.inc

@@ -45,7 +45,7 @@ subroutine da_transform_vtoy(cv_size, be, ep, cv, iv, vp, vv, xbx, &
 #ifdef VAR4D
       call da_transform_vtox(grid, be%cv%size_jb, xbx, be, ep, cv(1:be%cv%size_jb), vv, vp)
       if (be % ne > 0 .and. alphacv_method == alphacv_method_xa) then
-         call da_transform_vpatox(grid, be, ep, vp)
+         call da_transform_vpatox(grid, be%ne, ep, vp)
          call da_add_xa(grid%xa, grid%xa_ens) !grid%xa = grid%xa + xa_ens
       end if
       call domain_clock_get( grid, start_timestr=timestr )
@@ -195,7 +195,7 @@ subroutine da_transform_vtoy(cv_size, be, ep, cv, iv, vp, vv, xbx, &
 
       do nobwin = num_subtwindow, 1 , -1
          if (be % ne > 0 .and. alphacv_method == alphacv_method_xa) then
-            call da_transform_vpatox(grid, be, ep, vp, nobwin)
+            call da_transform_vpatox(grid, be%ne, ep, vp, nobwin)
             if ( num_subtwindow > 1 ) then
                !reset grid%xa = grid%xa_static
                call da_copy_xa(grid%xa, grid%xa_static)

var/da/da_minimisation/da_transform_vtoy_adj.inc

@@ -204,7 +204,7 @@ subroutine da_transform_vtoy_adj(cv_size, be, ep, cv, iv, vp, vv, xbx, y, &
       call da_transfer_xatowrftl_adj(grid, config_flags, 'gr01')
       !as of V3.9, da_zero_vp_type(vp) is commented out in the da_transform_vtox_adj,
       call da_zero_vp_type (vp) !as of V3.9, zero out vp here
-      call da_transform_vpatox_adj(grid, be, ep, vp)
+      call da_transform_vpatox_adj(grid, be%ne, ep, vp)
       call da_transform_vtox_adj(grid, be%cv%size_jb, xbx, be, ep, vp, vv, cv(1:be%cv%size_jb))
 
       if (num_fgat_time > 1 .and. use_rainobs) then
@@ -250,7 +250,7 @@ subroutine da_transform_vtoy_adj(cv_size, be, ep, cv, iv, vp, vv, xbx, y, &
 #endif
          call da_transform_xtoxa_adj(grid)
          call da_zero_vp_type (vp)
-         call da_transform_vpatox_adj(grid, be, ep, vp, nobwin)
+         call da_transform_vpatox_adj(grid, be%ne, ep, vp, nobwin)
          !as of V3.9, the calculation of vp%alpha is taken out of vptox_adj.
          !da_zero_vp_type(vp) is commented out in the da_transform_vtox_adj,
          !in order to pass in vp%alpha

var/da/da_setup_structures/da_setup_be_regional.inc

@@ -180,6 +180,12 @@ subroutine da_setup_be_regional(xb, be, grid)
 !===========for interpolating CV5==============================================================
 
    if (trace_use) call da_trace_entry("da_setup_be_regional")
+
+   ix = xb % mix
+   jy = xb % mjy
+   kz = xb % mkz
+
+if ( jb_factor > 0.0 ) then
    if( (cv_options == 5) .or. (cv_options == 6) .or. (cv_options == 7) ) then
       write (unit=message(1),fmt='(A,I3)') &
          'Set up background errors for regional application for cv_options = ',cv_options
@@ -204,10 +210,6 @@ subroutine da_setup_be_regional(xb, be, grid)
       call da_error(__FILE__,__LINE__,message(1:1))
    end if
 
-   ix = xb % mix
-   jy = xb % mjy
-   kz = xb % mkz
-
    be_rf_unit    = unit_end + 1
    be_print_unit = unit_end + 2
 
@@ -3617,6 +3619,23 @@ subroutine da_setup_be_regional(xb, be, grid)
    deallocate (bin)
    deallocate (bin2d)
 
+else ! jb_factor <= 0.0
+   ni = ix
+   nj = jy
+   nk = kz
+   be %  v1 % mz = 0
+   be %  v2 % mz = 0
+   be %  v3 % mz = 0
+   be %  v4 % mz = 0
+   be %  v5 % mz = 0
+   be %  v6 % mz = 0
+   be %  v7 % mz = 0
+   be %  v8 % mz = 0
+   be %  v9 % mz = 0
+   be % v10 % mz = 0
+   be % v11 % mz = 0
+end if ! jb_factor > 0.0
+
    if (be % ne > 0) then
       be % alpha % name = 'alpha'
       allocate (alpha_val(1:nk)) ! Not using jy dimension yet.

var/da/da_vtox_transforms/da_transform_vpatox.inc

@@ -1,20 +1,20 @@
-subroutine da_transform_vpatox(grid, be, ep, vp, nobwin) 
+subroutine da_transform_vpatox(grid, nens, ep, vp, nobwin)
 
    !----------------------------------------------------------------------
    ! Purpose:  Does transform of alpha control variable (vp%alpha) to x
    !----------------------------------------------------------------------
 
    implicit none
 
-   type(be_type),              intent(in)    :: be     ! background error structure.
    type(ep_type),              intent(in)    :: ep     ! Ensemble perturbation structure.
    type(vp_type),              intent(inout) :: vp     ! Grdipt/level CV.
    type(domain),               intent(inout) :: grid
+   integer, intent(in)                       :: nens   ! number of ensembles
    integer, intent(in), optional             :: nobwin
 
    integer :: iobwin
 
-   if (be%ne <= 0 .or. alphacv_method /= alphacv_method_xa) return
+   if (nens <= 0 .or. alphacv_method /= alphacv_method_xa) return
 
    if (trace_use) call da_trace_entry("da_transform_vpatox")
 
@@ -23,7 +23,7 @@ subroutine da_transform_vpatox(grid, be, ep, vp, nobwin)
       iobwin = nobwin
    end if
 
-   if (be % ne > 0 .and. alphacv_method == alphacv_method_xa) then
+   if (nens > 0 .and. alphacv_method == alphacv_method_xa) then
       !transform xa_ens from vp%alpha
       !xa_ens holds the increment from ensemble perturbations for certain fgat time
       grid%xa_ens%u = 0.0
@@ -37,9 +37,9 @@ subroutine da_transform_vpatox(grid, be, ep, vp, nobwin)
       grid%xa_ens%qsn = 0.0
       grid%xa_ens%qgr = 0.0
       if ( anal_type_hybrid_dual_res ) then
-         call da_calc_flow_dependence_xa_dual_res(grid, grid%xa_ens, be%ne, ep, vp, iobwin)
+         call da_calc_flow_dependence_xa_dual_res(grid, grid%xa_ens, nens, ep, vp, iobwin)
       else
-         call da_calc_flow_dependence_xa(grid%xa_ens, be%ne, ep, vp, iobwin)
+         call da_calc_flow_dependence_xa(grid%xa_ens, nens, ep, vp, iobwin)
       end if
    end if