wrfv2_fire/var/da/da_setup_structures/da_setup_be_regional.inc

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