var/da/da_setup_structures/da_setup_be_nmm_regional.inc

   1 subroutine da_setup_be_nmm_regional(xb, be)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: Define and allocate components of background errors
   5    !
   6    ! Updates:
   7    !
   8    !       Implementation of multi-variate BE
   9    !       Syed RH Rizvi,  MMM/NESL/NCAR,  Date: 02/01/2010
  10    !---------------------------------------------------------------------------
  11
  12    implicit none
  13
  14    type (xb_type), intent(in)    :: xb                    ! First guess structure.
  15    type (be_type), intent(inout) :: be                    ! Back. errors structure.
  16
  17    integer                     :: i, j, k, m       ! Loop counters.
  18    integer, allocatable:: bin(:,:,:)         ! Bin assigned to each 3D point
  19    integer, allocatable:: bin2d(:,:)         ! Bin assigned to each 2D point
  20    integer             :: bin_type           ! Type of bin to average over.
  21    integer             :: num_bins           ! Number of bins (3D fields).
  22    integer             :: num_bins2d         ! Number of bins (3D fields).
  23    real    :: lat_min, lat_max, binwidth_lat ! Used if bin_type = 2 (degrees)..
  24    real    :: hgt_min, hgt_max, binwidth_hgt ! Used if bin_type = 2 (m). .
  25
  26    real*8, allocatable         :: be1_eval_loc(:,:)     ! Temp arrays.
  27    real*8, allocatable         :: be2_eval_loc(:,:)     ! Temp arrays.
  28    real*8, allocatable         :: be3_eval_loc(:,:)     ! Temp arrays.
  29    real*8, allocatable         :: be4_eval_loc(:,:)     ! Temp arrays.
  30    real*8, allocatable         :: be5_eval_loc(:,:)     ! Temp arrays.
  31
  32    real*8, allocatable         :: be1_eval_glo(:)       ! Global Eigenvalues.
  33    real*8, allocatable         :: be2_eval_glo(:)       ! Global Eigenvalues.
  34    real*8, allocatable         :: be3_eval_glo(:)       ! Global Eigenvalues.
  35    real*8, allocatable         :: be4_eval_glo(:)       ! Global Eigenvalues.
  36    real*8, allocatable         :: be5_eval_glo(:)       ! Global Eigenvalues.
  37    real*8, allocatable         :: alpha_val(:)          ! Global Eigenvalues.
  38
  39    real*8, allocatable         :: be1_evec_loc(:,:,:)   ! Local Eigenvectors.
  40    real*8, allocatable         :: be2_evec_loc(:,:,:)   ! Local Eigenvectors.
  41    real*8, allocatable         :: be3_evec_loc(:,:,:)   ! Local Eigenvectors.
  42    real*8, allocatable         :: be4_evec_loc(:,:,:)   ! Local Eigenvectors.
  43    real*8, allocatable         :: be5_evec_loc(:,:,:)   ! Local Eigenvectors.
  44
  45    real*8, allocatable         :: be1_evec_glo(:,:)     ! Global Eigenvectors.
  46    real*8, allocatable         :: be2_evec_glo(:,:)     ! Global Eigenvectors.
  47    real*8, allocatable         :: be3_evec_glo(:,:)     ! Global Eigenvectors.
  48    real*8, allocatable         :: be4_evec_glo(:,:)     ! Global Eigenvectors.
  49    real*8, allocatable         :: be5_evec_glo(:,:)     ! Global Eigenvectors.
  50    real, allocatable           :: alpha_evec(:,:)       ! Global Eigenvectors.
  51
  52    real*8, allocatable         :: be1_rf_lengthscale(:) ! RF lengthscale.
  53    real*8, allocatable         :: be2_rf_lengthscale(:) ! RF lengthscale.
  54    real*8, allocatable         :: be3_rf_lengthscale(:) ! RF lengthscale.
  55    real*8, allocatable         :: be4_rf_lengthscale(:) ! RF lengthscale.
  56    real*8, allocatable         :: be5_rf_lengthscale(:)
  57    real*8, allocatable         :: alpha_rf_lengthscale(:)
  58    real*8, allocatable         :: alpha_rf_scale_factor(:)
  59
  60    real, allocatable           :: evec_loc(:,:,:)        ! Latitudinally varying eigenvectors.
  61    real, allocatable           :: eval_loc(:,:)          ! Latitudinally varying eigenvalues.
  62
  63    character*10                :: variable
  64    character*80                :: var80c
  65    integer                     :: ni, nj, nk, nk_2d, b
  66    integer                     :: ix, jy, kz, mz
  67    real, allocatable   :: regcoeff_psi_chi(:)        ! psi/chi    regression cooefficient.
  68    real, allocatable   :: regcoeff_psi_t(:,:,:)      ! psi/t   regression cooefficient.
  69    real, allocatable   :: regcoeff_psi_ps(:,:)       ! psi/ps     regression cooefficient.
  70    real, allocatable   :: regcoeff_psi_rh(:,:,:)     ! psi/rh     regression cooefficient.
  71    real, allocatable   :: regcoeff_chi_u_t(:,:,:)    ! chi_u/t regression coefficient
  72    real, allocatable   :: regcoeff_chi_u_ps(:,:)     ! chi_u/ps   regression coefficient
  73    real, allocatable   :: regcoeff_chi_u_rh(:,:,:)   ! chi_u/rh   regression coefficient
  74    real, allocatable   :: regcoeff_t_u_rh(:,:,:)     ! t_u/rh  regression coefficient
  75    real, allocatable   :: regcoeff_ps_u_rh(:,:)      ! ps_u/rh    regression coefficient
  76
  77    integer                     :: be_unit, ier, be_rf_unit, be_print_unit, it
  78
  79    if (trace_use) call da_trace_entry("da_setup_be_nmm_regional")
  80
  81    if( (cv_options == 5) .or. (cv_options == 6) .or. (cv_options == 7) ) then
  82       write (unit=message(1),fmt='(A,i3)') &
  83           'Set up background errors for regional application for cv_option= ',cv_options
  84       call da_message(message(1:1))
  85    else
  86       write(unit=message(1),fmt='(A,I4)') &
  87           'This subroutine is for cv_options = 5, 6, or 7; you have selected cv_options = ', cv_options
  88       call da_error(__FILE__,__LINE__,message(1:1))
  89    end if
  90
  91    ix = xb % mix
  92    jy = xb % mjy
  93    kz = xb % mkz
  94
  95    be_rf_unit    = unit_end + 1
  96    be_print_unit = unit_end + 2
  97
  98    call da_get_unit(be_unit)
  99    open(unit=be_unit,file="be.dat", status="old",form="unformatted")
 100
 101    rewind (be_unit)
 102    read (be_unit, iostat=ier) ni, nj, nk
 103    if (ier /= 0) then
 104       write(unit=message(1),fmt='(A,I4,A,I4)') 'cv_options:', cv_options,' Reading error in unit=',be_unit
 105       call da_error(__FILE__,__LINE__,message(1:1))
 106    else
 107       if ( nk /= kz ) then
 108          call da_error(__FILE__,__LINE__,  &
 109             (/"Vertical levels in fg and be.dat do not match."/))
 110       end if
 111    endif
 112
 113    allocate (bin(1:ni,1:nj,1:nk))
 114    allocate (bin2d(1:ni,1:nj))
 115
 116    read (be_unit)bin_type
 117    read (be_unit)lat_min, lat_max, binwidth_lat
 118    read (be_unit)hgt_min, hgt_max, binwidth_hgt
 119    read (be_unit)num_bins, num_bins2d
 120    read (be_unit)bin(1:ni,1:nj,1:nk)
 121    read (be_unit)bin2d(1:ni,1:nj)
 122
 123    ! 1.1 Read in regression coefficients
 124    allocate  (regcoeff_psi_chi(1:num_bins))
 125    allocate  (regcoeff_psi_t(1:nk,1:nk,1:num_bins2d))
 126    allocate  (regcoeff_psi_ps(1:nk,1:num_bins2d))
 127    allocate  (regcoeff_psi_rh(1:nk,1:nk,1:num_bins2d))
 128    allocate  (regcoeff_chi_u_t(1:nk,1:nk,1:num_bins2d))
 129    allocate  (regcoeff_chi_u_ps(1:nk,1:num_bins2d))
 130    allocate  (regcoeff_chi_u_rh(1:nk,1:nk,1:num_bins2d))
 131    allocate  (regcoeff_t_u_rh(1:nk,1:nk,1:num_bins2d))
 132    allocate  (regcoeff_ps_u_rh(1:nk,1:num_bins2d))
 133
 134    regcoeff_psi_chi    = 0.
 135    regcoeff_psi_t      = 0.
 136    regcoeff_psi_ps     = 0.
 137    regcoeff_psi_rh     = 0.
 138    regcoeff_chi_u_t    = 0.
 139    regcoeff_chi_u_ps   = 0.
 140    regcoeff_chi_u_rh   = 0.
 141    regcoeff_t_u_rh     = 0.
 142    regcoeff_ps_u_rh    = 0.
 143
 144    if ( (cv_options == 5) .or. (cv_options == 7) ) then
 145    read (be_unit) regcoeff_psi_chi
 146    read (be_unit) regcoeff_psi_ps
 147    read (be_unit) regcoeff_psi_t
 148    else
 149    do i = 1 , 9
 150    read (be_unit) var80c
 151    select case( trim(adjustl(var80c)) )
 152
 153    case ('regcoeff_psi_chi')
 154    read (be_unit) regcoeff_psi_chi
 155    case ('regcoeff_psi_t')
 156    read (be_unit) regcoeff_psi_t
 157    case ('regcoeff_psi_ps')
 158    read (be_unit) regcoeff_psi_ps
 159    case ('regcoeff_psi_rh')
 160    read (be_unit) regcoeff_psi_rh
 161    case ('regcoeff_chi_u_t')
 162    read (be_unit) regcoeff_chi_u_t
 163    case ('regcoeff_chi_u_ps')
 164    read (be_unit) regcoeff_chi_u_ps
 165    case ('regcoeff_chi_u_rh')
 166    read (be_unit) regcoeff_chi_u_rh
 167    case ('regcoeff_t_u_rh')
 168    read (be_unit) regcoeff_t_u_rh
 169    case ('regcoeff_ps_u_rh')
 170    read (be_unit) regcoeff_ps_u_rh
 171    case default;
 172       message(1)=' Read problem in regression coefficients in BE file '
 173       write (unit=message(2),fmt='(A,A)') ' Trying to read regression coefficients for variable: ',trim(adjustl(var80c))
 174       call da_error(__FILE__,__LINE__,message(1:2))
 175    end select
 176    end do
 177    end if
 178
 179    ! 1.2 Fill regression coeff. array
 180    allocate (be%reg_psi_chi (1:jy,1:nk))
 181    allocate (be%reg_psi_t   (1:jy,1:nk,1:nk))
 182    allocate (be%reg_psi_ps  (1:jy,1:nk))
 183    allocate (be%reg_psi_rh  (1:jy,1:nk,1:nk))
 184    allocate (be%reg_chi_u_t (1:jy,1:nk,1:nk))
 185    allocate (be%reg_chi_u_ps(1:jy,1:nk))
 186    allocate (be%reg_chi_u_rh(1:jy,1:nk,1:nk))
 187    allocate (be%reg_t_u_rh  (1:jy,1:nk,1:nk))
 188    allocate (be%reg_ps_u_rh (1:jy,1:nk))
 189
 190    be%reg_psi_chi = 0.
 191    be%reg_psi_t   = 0.
 192    be%reg_psi_ps  = 0.
 193    be%reg_psi_rh  = 0.
 194    be%reg_chi_u_t = 0.
 195    be%reg_chi_u_ps= 0.
 196    be%reg_chi_u_rh= 0.
 197    be%reg_t_u_rh  = 0.
 198
 199    do k=1,nk
 200       do j =1, jy
 201          b = bin(1,1,k)
 202          be%reg_psi_chi(j,k) = psi_chi_factor * regcoeff_psi_chi(b)
 203       end do
 204    end do
 205
 206    do j=1,jy
 207       b = bin2d(1,1)
 208       do k=1,nk
 209          be%reg_psi_ps(j,k)   = psi_ps_factor   * regcoeff_psi_ps(k,b)
 210          be%reg_ps_u_rh(j,k)  = ps_u_rh_factor  * regcoeff_ps_u_rh(k,b)
 211          be%reg_chi_u_ps(j,k) = chi_u_ps_factor * regcoeff_chi_u_ps(k,b)
 212       end do
 213    end do
 214
 215    do j=1,jy
 216       b = bin2d(1,1)
 217       do i=1,nk
 218          do k=1,nk
 219             be%reg_psi_t(j,i,k)   = psi_t_factor   *  regcoeff_psi_t(i,k,b)
 220             be%reg_psi_rh(j,i,k)  = psi_rh_factor  *  regcoeff_psi_rh(i,k,b)
 221             be%reg_chi_u_t(j,i,k) = chi_u_t_factor *  regcoeff_chi_u_t(i,k,b)
 222             be%reg_chi_u_rh(j,i,k)= chi_u_rh_factor*  regcoeff_chi_u_rh(i,k,b)
 223             be%reg_t_u_rh(j,i,k)  = t_u_rh_factor  *  regcoeff_t_u_rh(i,k,b)
 224          end do
 225       end do
 226    end do
 227
 228    deallocate (regcoeff_psi_chi)
 229    deallocate (regcoeff_psi_t)
 230    deallocate (regcoeff_psi_ps)
 231    deallocate (regcoeff_psi_rh)
 232    deallocate (regcoeff_chi_u_t)
 233    deallocate (regcoeff_chi_u_ps)
 234    deallocate (regcoeff_chi_u_rh)
 235    deallocate (regcoeff_t_u_rh)
 236    deallocate (regcoeff_ps_u_rh)
 237
 238    ! 1.3 Domain_averaged regression coefficients
 239
 240   if (.not.lat_stats_option) then
 241       write (unit=message(4), fmt='(3x, a)') &
 242          'Using the averaged regression coefficients for unbalanced part'
 243
 244       do k=1,nk
 245          be%reg_psi_ps  (1:num_bins2d,k) = sum(be%reg_psi_ps  (1:num_bins2d,k))/dble(num_bins2d)
 246          be%reg_ps_u_rh (1:num_bins2d,k) = sum(be%reg_ps_u_rh (1:num_bins2d,k))/dble(num_bins2d)
 247          be%reg_chi_u_ps(1:num_bins2d,k) = sum(be%reg_chi_u_ps(1:num_bins2d,k))/dble(num_bins2d)
 248       end do
 249
 250       do m=1,nk
 251          do k=1,nk
 252
 253            be%reg_psi_t  (1:num_bins2d,k,m)= sum(be%reg_psi_t(1:num_bins2d,k,m))/dble(num_bins2d)
 254            be%reg_psi_rh (1:num_bins2d,k,m)= sum(be%reg_psi_rh  (1:num_bins2d,k,m))/dble(num_bins2d)
 255            be%reg_chi_u_t(1:num_bins2d,k,m)= sum(be%reg_chi_u_t(1:num_bins2d,k,m))/dble(num_bins2d)
 256            be%reg_chi_u_rh(1:num_bins2d,k,m)= sum(be%reg_chi_u_rh  (1:num_bins2d,k,m))/dble(num_bins2d)
 257            be%reg_t_u_rh (1:num_bins2d,k,m)= sum(be%reg_t_u_rh (1:num_bins2d,k,m))/dble(num_bins2d)
 258
 259          end do
 260       end do
 261    else
 262       write (unit=message(4), fmt='(3x, a)') &
 263          'Using the latitude-dependent regression coefficients for unbalanced part'
 264    end if
 265
 266    call da_message(message(1:4))
 267
 268    ! 2.0 Load the eigenvector and eigenvalue
 269
 270    allocate (be1_eval_loc (1:jy,1:kz))
 271    allocate (be2_eval_loc (1:jy,1:kz))
 272    allocate (be3_eval_loc (1:jy,1:kz))
 273    allocate (be4_eval_loc (1:jy,1:kz))
 274    allocate (be5_eval_loc (1:jy,1:1))
 275
 276    if (vert_corr == vert_corr_2) then
 277
 278       allocate (be1_eval_glo(1:kz))
 279       allocate (be2_eval_glo(1:kz))
 280       allocate (be3_eval_glo(1:kz))
 281       allocate (be4_eval_glo(1:kz))
 282       allocate (be5_eval_glo(1:1))
 283
 284       allocate (be1_evec_loc(1:jy,1:kz,1:kz))
 285       allocate (be2_evec_loc(1:jy,1:kz,1:kz))
 286       allocate (be3_evec_loc(1:jy,1:kz,1:kz))
 287       allocate (be4_evec_loc(1:jy,1:kz,1:kz))
 288       allocate (be5_evec_loc(1:jy,1: 1,1: 1))
 289
 290       allocate (be1_evec_glo(1:kz,1:kz))
 291       allocate (be2_evec_glo(1:kz,1:kz))
 292       allocate (be3_evec_glo(1:kz,1:kz))
 293       allocate (be4_evec_glo(1:kz,1:kz))
 294       allocate (be5_evec_glo(1:1,1:1))
 295    end if
 296
 297    ! 2.2 Read in the eigenvector and eigenvalue
 298
 299    do i = 1 , 4
 300     read (be_unit) variable
 301    select case( trim(adjustl(variable)) )
 302    case ('psi')
 303    read (be_unit) nk, num_bins2d
 304    read (be_unit)  be1_evec_glo
 305    read (be_unit)  be1_eval_glo
 306    if( i == 1) then
 307    allocate (evec_loc(1:nk,1:nk,1:num_bins2d))
 308    allocate (eval_loc(1:nk,     1:num_bins2d))
 309    end if
 310    read (be_unit)  evec_loc
 311    read (be_unit)  eval_loc
 312    do j=1,jy
 313       b = bin2d(1,1)
 314       be1_evec_loc(j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
 315       be1_eval_loc(j,1:nk  ) = eval_loc(1:nk,b)
 316    end do
 317    be % v1 % name = trim(adjustl(variable))
 318
 319    case ('chi_u')
 320    read (be_unit) nk, num_bins2d
 321    read (be_unit)  be2_evec_glo
 322    read (be_unit)  be2_eval_glo
 323    read (be_unit)  evec_loc
 324    read (be_unit)  eval_loc
 325    do j=1,jy
 326       b = bin2d(1,1)
 327       be2_evec_loc(j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
 328       be2_eval_loc(j,1:nk  ) = eval_loc(1:nk,b)
 329    end do
 330    be % v2 % name = trim(adjustl(variable))
 331
 332    case ('t_u')
 333    read (be_unit) nk, num_bins2d
 334    read (be_unit)  be3_evec_glo
 335    read (be_unit)  be3_eval_glo
 336    read (be_unit)  evec_loc
 337    read (be_unit)  eval_loc
 338    do j=1,jy
 339       b = bin2d(1,1)
 340       be3_evec_loc(j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
 341       be3_eval_loc(j,1:nk  ) = eval_loc(1:nk,b)
 342    end do
 343    be % v3 % name = trim(adjustl(variable))
 344
 345    case ('rh_u' , 'rh' )
 346    read (be_unit) nk, num_bins2d
 347    read (be_unit)  be4_evec_glo
 348    read (be_unit)  be4_eval_glo
 349    read (be_unit)  evec_loc
 350    read (be_unit)  eval_loc
 351    do j=1,jy
 352       b = bin2d(1,1)
 353       be4_evec_loc(j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
 354       be4_eval_loc(j,1:nk  ) = eval_loc(1:nk,b)
 355    end do
 356    be % v4 % name = trim(adjustl(variable))
 357    case default;
 358       message(1)=' Read problem in eigen vaectors/values in BE file '
 359       write (unit=message(2),fmt='(A,A)') ' Trying to read Eigenvectors for variable: ',trim(adjustl(variable))
 360       call da_error(__FILE__,__LINE__,message(1:2))
 361    end select
 362    end do
 363
 364    deallocate (evec_loc)
 365    deallocate (eval_loc)
 366
 367    ! 2.2.5 Control variable ps_u
 368    read (be_unit) variable
 369    read (be_unit) nk_2d, num_bins2d
 370    allocate (evec_loc(1:nk_2d,1:nk_2d,1:num_bins2d))
 371    allocate (eval_loc(1:nk_2d,        1:num_bins2d))
 372    read (be_unit)  be5_evec_glo
 373    read (be_unit)  be5_eval_glo
 374    read (be_unit)  evec_loc
 375    read (be_unit)  eval_loc
 376    if( trim(adjustl(variable)) /= 'ps_u' ) then
 377       message(1)=' Variable mismatch while transfering eigen values from BE file '
 378       write (unit=message(2),fmt='(A,A)') ' Expected ps_u but got ',trim(adjustl(variable))
 379       call da_error(__FILE__,__LINE__,message(1:2))
 380    end if
 381    be % v5 % name = trim(adjustl(variable))
 382    do j=1,jy
 383       b = bin2d(1,1)
 384       be5_evec_loc(j,1:1,1:1) = evec_loc(1:1,1:1,b)
 385       be5_eval_loc(j,1:1 ) = eval_loc(1:1,b)
 386    end do
 387
 388    deallocate (evec_loc)
 389    deallocate (eval_loc)
 390
 391 !
 392    if(use_radarobs .and. use_radar_rf .or. use_rad .and. crtm_cloud) then
 393       if ( cloud_cv_options == 1 ) be % v4 % name = 'qt   '
 394    end if
 395
 396    write (unit=message(2),fmt='(3x,A,7A)') 'WRF-Var dry control variables are:', &
 397       trim(be % v1 % name), ', ', trim(be % v2 % name), ', ', &
 398       trim(be % v3 % name), ' and ', trim(be % v5 % name)
 399
 400    write (unit=message(3),fmt='(3x,A,A)') &
 401       'Humidity control variable is ', trim(be % v4 % name)
 402 !
 403    ! 3.0 Check and get the truncated number of the vertical modes
 404    ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 405
 406    if (vert_corr == vert_corr_2) then
 407
 408       ! 3.1 Perform checks on eigenvectors:
 409
 410       if (test_statistics) then
 411          call da_check_eof_decomposition(be1_eval_glo(:), be1_evec_glo(:,:), be % v1 % name)
 412          call da_check_eof_decomposition(be2_eval_glo(:), be2_evec_glo(:,:), be % v2 % name)
 413          call da_check_eof_decomposition(be3_eval_glo(:), be3_evec_glo(:,:), be % v3 % name)
 414          call da_check_eof_decomposition(be4_eval_glo(:), be4_evec_glo(:,:), be % v4 % name)
 415       end if
 416
 417       ! 3.2 Truncate in vertical:
 418
 419       call da_get_vertical_truncation(max_vert_var1, be1_eval_glo(:), be % v1)
 420       call da_get_vertical_truncation(max_vert_var2, be2_eval_glo(:), be % v2)
 421       call da_get_vertical_truncation(max_vert_var3, be3_eval_glo(:), be % v3)
 422       call da_get_vertical_truncation(max_vert_var4, be4_eval_glo(:), be % v4)
 423
 424       if (max_vert_var5 == 0.0) then
 425          be % v5 % mz = 0
 426       else
 427          be % v5 % mz = 1
 428       end if
 429
 430       write (unit=stdout,fmt=*) ' '
 431
 432    else
 433
 434       ! 3.3 no truncated
 435
 436       be % v1 % mz = xb % mkz
 437       be % v2 % mz = xb % mkz
 438       be % v3 % mz = xb % mkz
 439       be % v4 % mz = xb % mkz
 440       be % v5 % mz = xb % mkz
 441
 442    end if
 443
 444    ! 4.0 Initialise control variable space components of header:
 445    ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 446
 447    ! 4.1 Compute the size of the control variables
 448
 449    be % mix = ix
 450    be % mjy = jy
 451
 452    ! 4.2 Transfer errors to error structure:
 453
 454    call da_allocate_background_errors(jy, kz, be1_eval_glo, be1_evec_glo, be1_eval_loc, &
 455                                        be1_evec_loc, be % v1)
 456    call da_allocate_background_errors(jy, kz, be2_eval_glo, be2_evec_glo, be2_eval_loc, &
 457                                        be2_evec_loc, be % v2)
 458    call da_allocate_background_errors(jy, kz, be3_eval_glo, be3_evec_glo, be3_eval_loc, &
 459                                        be3_evec_loc, be % v3)
 460    call da_allocate_background_errors(jy, kz, be4_eval_glo, be4_evec_glo, be4_eval_loc, &
 461                                        be4_evec_loc, be % v4)
 462
 463    ! 4.2.1 transfer the ps_u variance to be % v5:
 464
 465    call da_allocate_background_errors(jy,  1, be5_eval_glo, be5_evec_glo, be5_eval_loc, &
 466                                        be5_evec_loc, be % v5)
 467
 468    if (print_detail_be) then
 469       write (unit=stdout,fmt='(3x,a,i10)') "b) Finished eigenvector processing!"
 470    end if
 471
 472    ! 5.0 Load the scale lengths
 473    ! ~~~~~~~~~~~~~~~~~~~~~~~~~~
 474
 475    ! 5.1 allocate the array for scale lengths
 476
 477    allocate (be1_rf_lengthscale(1:nk))
 478    allocate (be2_rf_lengthscale(1:nk))
 479    allocate (be3_rf_lengthscale(1:nk))
 480    allocate (be4_rf_lengthscale(1:nk))
 481    allocate (be5_rf_lengthscale(1:nk))
 482
 483    ! 5.2 read in the scale lengths
 484    do i = 1 , 5
 485     read (be_unit) variable
 486    select case( trim(adjustl(variable)) )
 487
 488    case ('psi')
 489      read(be_unit) be1_rf_lengthscale
 490    case ('chi_u')
 491      read(be_unit) be2_rf_lengthscale
 492    case ('t_u')
 493      read(be_unit) be3_rf_lengthscale
 494    case ('rh_u' , 'rh')
 495      read(be_unit) be4_rf_lengthscale
 496    case ('ps_u')
 497      read(be_unit) be5_rf_lengthscale(1:1)
 498    case default;
 499       message(1)=' Read problem in lengthscales in BE file '
 500       write (unit=message(2),fmt='(A,A)') ' Trying to read lengthscales for variable: ',trim(adjustl(variable))
 501       call da_error(__FILE__,__LINE__,message(1:2))
 502    end select
 503    end do
 504
 505
 506    ! 5.3 Convert the scale lengths in the real distance (meter)
 507
 508    be1_rf_lengthscale(1:nk) = be1_rf_lengthscale(1:nk) * xb%ds
 509    be2_rf_lengthscale(1:nk) = be2_rf_lengthscale(1:nk) * xb%ds
 510    be3_rf_lengthscale(1:nk) = be3_rf_lengthscale(1:nk) * xb%ds
 511    be4_rf_lengthscale(1:nk) = be4_rf_lengthscale(1:nk) * xb%ds
 512    be5_rf_lengthscale(1:1)  = be5_rf_lengthscale(1:1)  * xb%ds
 513
 514    ! 6.0 Perform checks on eigenvectors with be data structure:
 515    if (test_statistics) then
 516       call da_check_eof_decomposition(be%v1%val_g(:), be%v1%evec_g(:,:),&
 517                                      be%v1%name)
 518       call da_check_eof_decomposition(be%v2%val_g(:), be%v2%evec_g(:,:),&
 519                                      be%v2%name)
 520       call da_check_eof_decomposition(be%v3%val_g(:), be%v3%evec_g(:,:),&
 521                                      be%v3%name)
 522       call da_check_eof_decomposition(be%v4%val_g(:), be%v4%evec_g(:,:),&
 523                                      be%v4%name)
 524    end if
 525
 526    ! 6.2 Close the be unit
 527
 528    close(be_unit)
 529    call da_free_unit(be_unit)
 530
 531    ! 6.3 Keep the original be % v1, be % v2,...., and lengthscale in the first loop
 532    !     for the rescaling in the later loops:
 533
 534    it = 1
 535
 536    if (max_ext_its > 1) then
 537
 538        write(unit=message(1),fmt='(A,I4)') '>>> Save the variances and scale-lengths in outer-loop', it
 539        call da_message(message(1:1))
 540        write(be_rf_unit)  kz, jy, ix, be % v1 % mz, be % v2 % mz, be% v3 % mz, &
 541                                      be % v4 % mz, be % v5 % mz, xb % ds
 542        write(be_rf_unit) be % v1 % val, be % v2 % val, be% v3 % val, &
 543                                      be % v4 % val, be % v5 % val, &
 544              be1_rf_lengthscale, be2_rf_lengthscale, be3_rf_lengthscale, &
 545              be4_rf_lengthscale, be5_rf_lengthscale
 546
 547        if (print_detail_be ) then
 548        write(be_print_unit,'("it=",i2,2x,"kz=",i3,2x,"jy=",i4,2x,"ix=",i4,2x,"ds=",e12.5)') &
 549                                                it, kz, jy, ix, xb % ds
 550        write(be_print_unit,'("Original val and rf, and mz:",5i5)') &
 551                       be % v1 % mz, be % v2 % mz, be% v3 % mz, be % v4 % mz, be % v5 % mz
 552        write(be_print_unit,'("mz=",i3,2x,"be%v1%val:"/(10e12.5))') be%v1%mz, be%v1%val(1,:)
 553        write(be_print_unit,'("mz=",i3,2x,"be%v2%val:"/(10e12.5))') be%v2%mz, be%v2%val(1,:)
 554        write(be_print_unit,'("mz=",i3,2x,"be%v3%val:"/(10e12.5))') be%v3%mz, be%v3%val(1,:)
 555        write(be_print_unit,'("mz=",i3,2x,"be%v4%val:"/(10e12.5))') be%v4%mz, be%v4%val(1,:)
 556        write(be_print_unit,'("mz=",i3,2x,"be%v5%val:"/(10e12.5))') be%v5%mz, be%v5%val(1,:)
 557        write(be_print_unit,'(/"scale-length: kz=",i3)') kz
 558        do i = 1,kz
 559          if (i == 1) then
 560            write(be_print_unit,'(i3,2x,5e15.5)') i,be1_rf_lengthscale(i), &
 561                 be2_rf_lengthscale(i), be3_rf_lengthscale(i), be4_rf_lengthscale(i), &
 562                 be5_rf_lengthscale(i)
 563          else
 564            write(be_print_unit,'(i3,2x,5e15.5)') i,be1_rf_lengthscale(i), &
 565                 be2_rf_lengthscale(i), be3_rf_lengthscale(i), be4_rf_lengthscale(i)
 566        endif
 567        enddo
 568
 569        endif
 570
 571    endif
 572
 573    ! 7.0 Apply empirical and recursive filter rescaling factor:
 574    ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 575
 576    call da_rescale_background_errors(var_scaling1(1), len_scaling1(1), &
 577                                       xb % ds, be1_rf_lengthscale, be % v1)
 578    call da_rescale_background_errors(var_scaling2(1), len_scaling2(1), &
 579                                       xb % ds, be2_rf_lengthscale, be % v2)
 580    call da_rescale_background_errors(var_scaling3(1), len_scaling3(1), &
 581                                       xb % ds, be3_rf_lengthscale, be % v3)
 582    call da_rescale_background_errors(var_scaling4(1), len_scaling4(1), &
 583                                       xb % ds, be4_rf_lengthscale, be % v4)
 584    call da_rescale_background_errors(var_scaling5(1), len_scaling5(1), &
 585                                       xb % ds, be5_rf_lengthscale, be % v5)
 586    if (print_detail_be ) then
 587
 588        write(be_print_unit,'("it=",i2,2x,"kz=",i3,2x,"jy=",i4,2x,"ix=",i4,2x,"ds=",e12.5)') &
 589                                                       it, kz, jy, ix, xb % ds
 590        write(be_print_unit,'("Loop it=",i2," val and rf, and mz:",5i5)') &
 591                   it, be % v1 % mz, be % v2 % mz, be% v3 % mz, be % v4 % mz, be % v5 % mz
 592        write(be_print_unit,'("mz=",i3,2x,"be%v1%val:"/(10e12.5))') be%v1%mz, be%v1%val(1,:)
 593        write(be_print_unit,'("mz=",i3,2x,"be%v2%val:"/(10e12.5))') be%v2%mz, be%v2%val(1,:)
 594        write(be_print_unit,'("mz=",i3,2x,"be%v3%val:"/(10e12.5))') be%v3%mz, be%v3%val(1,:)
 595        write(be_print_unit,'("mz=",i3,2x,"be%v4%val:"/(10e12.5))') be%v4%mz, be%v4%val(1,:)
 596        write(be_print_unit,'("mz=",i3,2x,"be%v5%val:"/(10e12.5))') be%v5%mz, be%v5%val(1,:)
 597        write(be_print_unit,'(/"scale-length: kz=",i3)') kz
 598        do i = 1,kz
 599          if (i == 1) then
 600            write(be_print_unit,'(i3,2x,5e15.5)') i, be % v1 % rf_alpha(i), &
 601                 be % v2 % rf_alpha(i), be % v3 % rf_alpha(i), be % v4 % rf_alpha(i), &
 602                 be % v5 % rf_alpha(i)
 603          else
 604            write(be_print_unit,'(i3,2x,5e15.5)') i, be % v1 % rf_alpha(i), &
 605                 be % v2 % rf_alpha(i), be % v3 % rf_alpha(i), be % v4 % rf_alpha(i)
 606        endif
 607        enddo
 608
 609    endif
 610
 611    ! 8.0 deallocate input model state:
 612    ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 613
 614    deallocate (be1_eval_loc)
 615    deallocate (be2_eval_loc)
 616    deallocate (be3_eval_loc)
 617    deallocate (be4_eval_loc)
 618    deallocate (be5_eval_loc)
 619
 620    deallocate (be1_rf_lengthscale)
 621    deallocate (be2_rf_lengthscale)
 622    deallocate (be3_rf_lengthscale)
 623    deallocate (be4_rf_lengthscale)
 624    deallocate (be5_rf_lengthscale)
 625
 626    if (vert_corr == vert_corr_2) then
 627       deallocate (be1_eval_glo)
 628       deallocate (be2_eval_glo)
 629       deallocate (be3_eval_glo)
 630       deallocate (be4_eval_glo)
 631       deallocate (be5_eval_glo)
 632
 633       deallocate (be1_evec_loc)
 634       deallocate (be2_evec_loc)
 635       deallocate (be3_evec_loc)
 636       deallocate (be4_evec_loc)
 637       deallocate (be5_evec_loc)
 638
 639       deallocate (be1_evec_glo)
 640       deallocate (be2_evec_glo)
 641       deallocate (be3_evec_glo)
 642       deallocate (be4_evec_glo)
 643       deallocate (be5_evec_glo)
 644
 645    end if
 646
 647    deallocate (bin)
 648    deallocate (bin2d)
 649
 650    if (be % ne > 0) then
 651       be % alpha % name = 'alpha'
 652       allocate (alpha_val(1:kz)) ! Not using jy dimension yet.
 653       allocate (alpha_evec(1:kz,1:kz)) ! Not using jy dimension yet.
 654
 655       if ( alpha_vertloc_opt > 0 ) then ! Use vertical localization:
 656          call da_get_unit(be_unit)
 657          open(unit=be_unit,file='be.vertloc.dat', status='old', form='unformatted')
 658
 659          read (be_unit)nk
 660          read (be_unit)alpha_val(1:nk)
 661          read (be_unit)alpha_evec(1:nk,1:nk)
 662          close(be_unit)
 663
 664          call da_get_vertical_truncation(max_vert_var_alpha, alpha_val, be % alpha)
 665       else
 666          be % alpha % mz = 1 ! No vertical localization.
 667          alpha_val(1) = 1.0
 668          alpha_val(2:kz) = 0.0
 669          alpha_evec(:,:) = 1.0
 670       end if
 671       mz = be % alpha % mz
 672
 673 !     Alpha eigenvalues and eigenvectors:
 674       allocate (be % alpha % val(1:jy,1:mz)) ! Not using jy dimension but here for consistency.
 675       allocate (be % alpha % evec(1:jy,1:kz,1:mz))
 676       do m = 1, mz
 677          be % alpha % val(:,m) = sigma_alpha * alpha_val(m)
 678          do k = 1, nk
 679             be % alpha % evec(:,k,m) = alpha_evec(k,m)
 680          end do
 681       end do
 682
 683 !     Alpha RF lengthscales and variance scaling factors:
 684       allocate (alpha_rf_lengthscale(1:mz))
 685       allocate (be % alpha % rf_alpha(1:mz))
 686       allocate (alpha_rf_scale_factor(1:mz))
 687
 688       alpha_rf_lengthscale(1:mz) = 1000.0 * alpha_corr_scale / xb % ds ! Convert km to grid spacings.
 689
 690       call da_calculate_rf_factors( alpha_rf_lengthscale(:), be % alpha % rf_alpha(:), &
 691                                     alpha_rf_scale_factor(:) )
 692       do m = 1, mz
 693          be % alpha % val(:,m) = alpha_rf_scale_factor(m) * be % alpha % val(:,m)
 694       end do
 695
 696    else
 697       be % alpha % mz = 0
 698    end if
 699
 700    if (trace_use) call da_trace_exit("da_setup_be_nmm_regional")
 701
 702 end subroutine da_setup_be_nmm_regional
 703
 704