wrfv2_fire/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
   7    implicit none
   8
   9    type (xb_type), intent(in)  :: xb                    ! First guess structure.
  10    type (be_type), intent(out) :: 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
  33    real, allocatable           :: be1_evec_loc(:,:,:)   ! Local Eigenvectors.
  34    real, allocatable           :: be2_evec_loc(:,:,:)   ! Local Eigenvectors.
  35    real, allocatable           :: be3_evec_loc(:,:,:)   ! Local Eigenvectors.
  36    real, allocatable           :: be4_evec_loc(:,:,:)   ! Local Eigenvectors.
  37    real, allocatable           :: be5_evec_loc(:,:,:)   ! Local Eigenvectors.
  38
  39    real, allocatable           :: be1_evec_glo(:,:)     ! Global Eigenvectors.
  40    real, allocatable           :: be2_evec_glo(:,:)     ! Global Eigenvectors.
  41    real, allocatable           :: be3_evec_glo(:,:)     ! Global Eigenvectors.
  42    real, allocatable           :: be4_evec_glo(:,:)     ! Global Eigenvectors.
  43    real, allocatable           :: be5_evec_glo(:,:)     ! Global Eigenvectors.
  44
  45    real, allocatable           :: be1_rf_lengthscale(:) ! RF lengthscale.
  46    real, allocatable           :: be2_rf_lengthscale(:) ! RF lengthscale.
  47    real, allocatable           :: be3_rf_lengthscale(:) ! RF lengthscale.
  48    real, allocatable           :: be4_rf_lengthscale(:) ! RF lengthscale.
  49    real, allocatable           :: be5_rf_lengthscale(:)
  50    real, allocatable           :: alpha_rf_lengthscale(:)
  51
  52    real, allocatable           :: evec_loc(:,:,:)        ! Latitudinally varying eigenvectors.
  53    real, allocatable           :: eval_loc(:,:)          ! Latitudinally varying eigenvalues.
  54
  55    character*10                :: variable
  56    integer                     :: ni, nj, nk, nk_2d, b
  57    integer                     :: ix, jy, kz
  58    real, allocatable           :: regcoeff1(:)
  59    real, allocatable           :: regcoeff2(:,:)
  60    real, allocatable           :: regcoeff3(:,:,:)
  61    real                        :: avg,avg2,avg3
  62    integer                     :: be_unit
  63
  64    if (trace_use) call da_trace_entry("da_setup_be_nmm_regional")
  65
  66    be % v1 % name = 'psi  '           ! Streamfunction
  67    be % v2 % name = 'chi_u'           ! Uncorrelated velocity potential.
  68    be % v3 % name = 't_u'             ! Unbalanced temperature.
  69    be % v4 % name = 'q/qsg'
  70    be % v5 % name = 'psfc'            ! surface pressure
  71 !clq
  72    if(use_radarobs .and. use_radar_rf) then
  73    be % v4 % name = 'qt   '
  74    end if
  75 !clq end
  76
  77    write (unit=message(2),fmt='(3x,A,7A)') 'WRF-Var dry control variables are:', &
  78       trim(be % v1 % name), ', ', trim(be % v2 % name), ', ', &
  79       trim(be % v3 % name), ' and ', trim(be % v5 % name)
  80
  81    write (unit=message(3),fmt='(3x,A,A)') &
  82       'Humidity control variable is ', trim(be % v4 % name)
  83
  84    write (unit=message(3),fmt='(3x,A,A)') &
  85       'Humidity control variable is ', trim(be % v4 % name)
  86
  87    ix = xb % mix
  88    jy = xb % mjy
  89    kz = xb % mkz
  90    call da_get_unit(be_unit)
  91    open(unit=be_unit,file="be.dat", status="old",form="unformatted")
  92
  93    rewind (be_unit)
  94    read (be_unit) ni, nj, nk
  95
  96    allocate (bin(1:ni,1:nj,1:nk))
  97    allocate (bin2d(1:ni,1:nj))
  98
  99    read (be_unit)bin_type
 100    read (be_unit)lat_min, lat_max, binwidth_lat
 101    read (be_unit)hgt_min, hgt_max, binwidth_hgt
 102    read (be_unit)num_bins, num_bins2d
 103    read (be_unit)bin(1:ni,1:nj,1:nk)
 104    read (be_unit)bin2d(1:ni,1:nj)
 105
 106    ! 1.1 Read in regression coefficients
 107
 108
 109    allocate (regcoeff1(1:num_bins))
 110    allocate (regcoeff2(1:nk,1:num_bins2d))
 111    allocate (regcoeff3(1:nk,1:nk,1:num_bins2d))
 112
 113    read (be_unit) regcoeff1
 114    read (be_unit) regcoeff2
 115    read (be_unit) regcoeff3
 116
 117    ! 1.2 Fill regression coeff. array
 118
 119 !                                47 25 50
 120    allocate (be%reg_chi(1:jy,1:kz))
 121    allocate (be%reg_ps (1:jy,1:kz))
 122    allocate (be%reg_t  (1:jy,1:kz,1:kz))
 123
 124    do k=1,kz
 125       do j =1, jy
 126          b = bin(1,1,k)
 127          be%reg_chi(j,k) = regcoeff1(b)
 128       end do
 129    end do
 130    do j=1,jy
 131       b = bin2d(1,1)
 132       do k=1,kz
 133          be%reg_ps(j,k) = regcoeff2(k,b)
 134       end do
 135    end do
 136
 137    do j=1,jy
 138       b = bin2d(1,1)
 139       do i=1,kz
 140          do k=1,kz
 141             be%reg_t(j,i,k) = regcoeff3(i,k,b)
 142          end do
 143       end do
 144    end do
 145
 146    deallocate (regcoeff1)
 147    deallocate (regcoeff2)
 148    deallocate (regcoeff3)
 149
 150
 151    ! 1.3 Domain_averaged regression coefficients
 152
 153    if (.not.lat_stats_option) then
 154       write (unit=message(4), fmt='(3x, a)') &
 155          'Using the averaged regression coefficients for unbalanced part'
 156
 157       do k=1,nk
 158          avg= 0.0
 159          avg2=0.0
 160          do j=1,num_bins2d
 161             avg= avg + be%reg_ps (j,k)/float(num_bins2d)
 162             avg2= avg2 + be%reg_chi (j,k)/float(num_bins2d)
 163          end do
 164
 165          do j=1,num_bins2d
 166             be%reg_ps(j,k)=avg
 167             be%reg_chi (j,k)=avg2
 168          end do
 169       end do
 170
 171       do m=1,nk
 172          do k=1,nk
 173             avg3= 0.0
 174
 175             do j=1,num_bins2d
 176                avg3= avg3 + be%reg_t (j,k,m)/float(num_bins2d)
 177             end do
 178
 179             do j=1,num_bins2d
 180                be%reg_t(j,k,m)=avg3
 181             end do
 182          end do
 183       end do
 184    else
 185       write (unit=message(4), fmt='(3x, a)') &
 186          'Using the latitude-dependent regression coefficients for unbalanced part'
 187    end if
 188
 189    call da_message(message(1:4))
 190
 191    ! 2.0 Load the eigenvector and eigenvalue
 192    allocate (be1_eval_loc (1:jy,1:nk))
 193    allocate (be2_eval_loc (1:jy,1:nk))
 194    allocate (be3_eval_loc (1:jy,1:nk))
 195    allocate (be4_eval_loc (1:jy,1:nk))
 196    allocate (be5_eval_loc (1:jy,1:1))
 197
 198    if (vert_corr == vert_corr_2) then
 199
 200       allocate (be1_eval_glo(1:nk))
 201       allocate (be2_eval_glo(1:nk))
 202       allocate (be3_eval_glo(1:nk))
 203       allocate (be4_eval_glo(1:nk))
 204       allocate (be5_eval_glo(1:1))
 205
 206       allocate (be1_evec_loc(1:jy,1:nk,1:nk))
 207       allocate (be2_evec_loc(1:jy,1:nk,1:nk))
 208       allocate (be3_evec_loc(1:jy,1:nk,1:nk))
 209       allocate (be4_evec_loc(1:jy,1:nk,1:nk))
 210       allocate (be5_evec_loc(1:jy,1: 1,1: 1))
 211
 212       allocate (be1_evec_glo(1:nk,1:nk))
 213       allocate (be2_evec_glo(1:nk,1:nk))
 214       allocate (be3_evec_glo(1:nk,1:nk))
 215       allocate (be4_evec_glo(1:nk,1:nk))
 216       allocate (be5_evec_glo(1:1,1:1))
 217    end if
 218
 219    ! 2.2 Read in the eigenvector and eigenvalue
 220
 221    ! 2.2.1 Control variable 1 (psi)
 222
 223    read (be_unit) variable
 224    read (be_unit) nk, num_bins2d
 225    read (be_unit)  be1_evec_glo
 226    read (be_unit)  be1_eval_glo
 227    allocate (evec_loc(1:nk,1:nk,1:num_bins2d))
 228    allocate (eval_loc(1:nk,     1:num_bins2d))
 229    read (be_unit)  evec_loc
 230    read (be_unit)  eval_loc
 231    if (variable(1:3) /= 'psi') then
 232       message(1)=' Variable mismatch while transfering eigen values from BE file '
 233       write (unit=message(2),fmt='(A,A)') ' Expected psi but got ',variable
 234       call da_error(__FILE__,__LINE__,message(1:2))
 235    end if
 236    be % v1 % name = variable
 237    do j=1,jy
 238       b = bin2d(1,1)
 239       be1_evec_loc(j,1:kz,1:kz) = evec_loc(1:kz,1:kz,b)
 240       be1_eval_loc(j,1:kz  ) = eval_loc(1:kz,b)
 241    end do
 242
 243    ! 2.2.2 Control variable 2 (chi_u)
 244
 245    read (be_unit) variable
 246    read (be_unit) nk, num_bins2d
 247    read (be_unit)  be2_evec_glo
 248    read (be_unit)  be2_eval_glo
 249    read (be_unit)  evec_loc
 250    read (be_unit)  eval_loc
 251    if (variable(1:5) /= 'chi_u') then
 252       message(1)=' Variable mismatch while transfering eigen values from BE file '
 253       write (unit=message(2),fmt='(A,A)') ' Expected chi_u but got ',variable
 254       call da_error(__FILE__,__LINE__,message(1:2))
 255    end if
 256    be % v2 % name = variable
 257    do j=1,jy
 258       b = bin2d(1,1)
 259       be2_evec_loc(j,1:kz,1:kz) = evec_loc(1:kz,1:kz,b)
 260       be2_eval_loc(j,1:kz  ) = eval_loc(1:kz,b)
 261    end do
 262
 263    ! 2.2.3 Control variable 3 (t_u)
 264    read (be_unit) variable
 265    read (be_unit) nk, num_bins2d
 266    read (be_unit)  be3_evec_glo
 267    read (be_unit)  be3_eval_glo
 268    read (be_unit)  evec_loc
 269    read (be_unit)  eval_loc
 270    if (variable(1:3) /= 't_u') then
 271       message(1)=' Variable mismatch while transfering eigen values from BE file '
 272       write (unit=message(2),fmt='(A,A)') ' Expected t_u but got ',variable
 273       call da_error(__FILE__,__LINE__,message(1:2))
 274    end if
 275    be % v3 % name = variable
 276    do j=1,jy
 277       b = bin2d(1,1)
 278       be3_evec_loc(j,1:kz,1:kz) = evec_loc(1:kz,1:kz,b)
 279       be3_eval_loc(j,1:kz  ) = eval_loc(1:kz,b)
 280    end do
 281
 282    ! 2.2.4 Control variable 4 (q/qsg)
 283
 284    read (be_unit) variable
 285    read (be_unit) nk, num_bins2d
 286    read (be_unit)  be4_evec_glo
 287    read (be_unit)  be4_eval_glo
 288    read (be_unit)  evec_loc
 289    read (be_unit)  eval_loc
 290    if (variable(1:3) /= 'rh') then
 291       message(1)=' Variable mismatch while transfering eigen values from BE file '
 292       write (unit=message(2),fmt='(A,A)') ' Expected rh but got ',variable
 293       call da_error(__FILE__,__LINE__,message(1:2))
 294    end if
 295    be % v4 % name = variable
 296    do j=1,jy
 297       b = bin2d(1,1)
 298       be4_evec_loc(j,1:kz,1:kz) = evec_loc(1:kz,1:kz,b)
 299       be4_eval_loc(j,1:kz  ) = eval_loc(1:kz,b)
 300    end do
 301
 302    deallocate (evec_loc)
 303    deallocate (eval_loc)
 304
 305    ! 2.2.5 Control variable ps_u
 306
 307    read (be_unit) variable
 308    read (be_unit) nk_2d, num_bins2d
 309    allocate (evec_loc(1:nk_2d,1:nk_2d,1:num_bins2d))
 310    allocate (eval_loc(1:nk_2d,        1:num_bins2d))
 311    read (be_unit)  be5_evec_glo
 312    read (be_unit)  be5_eval_glo
 313    read (be_unit)  evec_loc
 314    read (be_unit)  eval_loc
 315    if (variable(1:4) /= 'ps_u') then
 316       message(1)=' Variable mismatch while transfering eigen values from BE file '
 317       write (unit=message(2),fmt='(A,A)') ' Expected ps_u but got ',variable
 318       call da_error(__FILE__,__LINE__,message(1:2))
 319    end if
 320    be % v5 % name = variable
 321    do j=1,jy
 322       b = bin2d(1,1)
 323       be5_evec_loc(j,1:1,1:1) = evec_loc(1:1,1:1,b)
 324       be5_eval_loc(j,1:1 ) = eval_loc(1:1,b)
 325    end do
 326
 327    deallocate (evec_loc)
 328    deallocate (eval_loc)
 329
 330    ! 3.0 Check and get the truncated number of the vertical modes
 331    ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 332
 333    if (vert_corr == vert_corr_2) then
 334
 335       ! 3.1 Perform checks on eigenvectors:
 336
 337       if (test_statistics) then
 338          call da_check_eof_decomposition(be1_eval_glo(:), be1_evec_glo(:,:), be % v1 % name)
 339          call da_check_eof_decomposition(be2_eval_glo(:), be2_evec_glo(:,:), be % v2 % name)
 340          call da_check_eof_decomposition(be3_eval_glo(:), be3_evec_glo(:,:), be % v3 % name)
 341          call da_check_eof_decomposition(be4_eval_glo(:), be4_evec_glo(:,:), be % v4 % name)
 342       end if
 343
 344       ! 3.2 Truncate in vertical:
 345
 346       call da_get_vertical_truncation(max_vert_var1, be1_eval_glo(:), be % v1)
 347       call da_get_vertical_truncation(max_vert_var2, be2_eval_glo(:), be % v2)
 348       call da_get_vertical_truncation(max_vert_var3, be3_eval_glo(:), be % v3)
 349       call da_get_vertical_truncation(max_vert_var4, be4_eval_glo(:), be % v4)
 350
 351       be % v5 % mz = 1
 352
 353       write (unit=stdout,fmt=*) ' '
 354
 355    else
 356
 357       ! 3.3 no truncated
 358
 359       be % v1 % mz = xb % mkz
 360       be % v2 % mz = xb % mkz
 361       be % v3 % mz = xb % mkz
 362       be % v4 % mz = xb % mkz
 363       be % v5 % mz = xb % mkz
 364
 365    end if
 366
 367    ! 4.0 Initialise control variable space components of header:
 368    ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 369
 370    ! 4.1 Compute the size of the control variables
 371
 372    be % mix = ix
 373    be % mjy = jy
 374
 375    ! 4.2 Transfer errors to error structure:
 376
 377    call da_allocate_background_errors(jy, kz, be1_eval_glo, be1_evec_glo, be1_eval_loc, &
 378                                        be1_evec_loc, be % v1)
 379    call da_allocate_background_errors(jy, kz, be2_eval_glo, be2_evec_glo, be2_eval_loc, &
 380                                        be2_evec_loc, be % v2)
 381    call da_allocate_background_errors(jy, kz, be3_eval_glo, be3_evec_glo, be3_eval_loc, &
 382                                        be3_evec_loc, be % v3)
 383    call da_allocate_background_errors(jy, kz, be4_eval_glo, be4_evec_glo, be4_eval_loc, &
 384                                        be4_evec_loc, be % v4)
 385    ! 4.2.1 transfer the ps_u variance to be % v5:
 386
 387    call da_allocate_background_errors(jy,  1, be5_eval_glo, be5_evec_glo, be5_eval_loc, &
 388                                        be5_evec_loc, be % v5)
 389
 390    if (print_detail_be) then
 391       write (unit=stdout,fmt='(3x,a,i10)') "b) Finished eigenvector processing!"
 392    end if
 393
 394    ! 5.0 Load the scale lengths
 395    ! ~~~~~~~~~~~~~~~~~~~~~~~~~~
 396
 397    ! 5.1 allocate the array for scale lengths
 398
 399    allocate (be1_rf_lengthscale(1:nk))
 400    allocate (be2_rf_lengthscale(1:nk))
 401    allocate (be3_rf_lengthscale(1:nk))
 402    allocate (be4_rf_lengthscale(1:nk))
 403    allocate (be5_rf_lengthscale(1:nk))
 404
 405    ! 5.2 read in the scale lengths
 406
 407    read (be_unit) variable
 408    read (be_unit) be1_rf_lengthscale
 409
 410    read (be_unit) variable
 411    read (be_unit) be2_rf_lengthscale
 412
 413    read (be_unit) variable
 414    read (be_unit) be3_rf_lengthscale
 415
 416    read (be_unit) variable
 417    read (be_unit) be4_rf_lengthscale
 418
 419    read (be_unit) variable
 420    read (be_unit) be5_rf_lengthscale(1:1)
 421    be%v5%name = variable
 422
 423    ! 5.3 Convert the scale lengths in the real distance (meter)
 424
 425    be1_rf_lengthscale(1:nk) = be1_rf_lengthscale(1:nk) * xb%ds
 426    be2_rf_lengthscale(1:nk) = be2_rf_lengthscale(1:nk) * xb%ds
 427    be3_rf_lengthscale(1:nk) = be3_rf_lengthscale(1:nk) * xb%ds
 428    be4_rf_lengthscale(1:nk) = be4_rf_lengthscale(1:nk) * xb%ds
 429    be5_rf_lengthscale(1:1)  = be5_rf_lengthscale(1:1)  * xb%ds
 430
 431    ! 6.0 Perform checks on eigenvectors with be data structure:
 432    ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 433
 434    if (test_statistics) then
 435       call da_check_eof_decomposition(be%v1%val_g(:), be%v1%evec_g(:,:),&
 436                                      be%v1%name)
 437       call da_check_eof_decomposition(be%v2%val_g(:), be%v2%evec_g(:,:),&
 438                                      be%v2%name)
 439       call da_check_eof_decomposition(be%v3%val_g(:), be%v3%evec_g(:,:),&
 440                                      be%v3%name)
 441       call da_check_eof_decomposition(be%v4%val_g(:), be%v4%evec_g(:,:),&
 442                                      be%v4%name)
 443    end if
 444
 445    ! 6.2 Close the be unit
 446
 447    close(be_unit)
 448    call da_free_unit(be_unit)
 449
 450    ! 7.0 Apply empirical and recursive filter rescaling factor:
 451    ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 452
 453    call da_rescale_background_errors(var_scaling1(1), len_scaling1(1), &
 454                                       xb % ds, be1_rf_lengthscale, be % v1)
 455    call da_rescale_background_errors(var_scaling2(1), len_scaling2(1), &
 456                                       xb % ds, be2_rf_lengthscale, be % v2)
 457    call da_rescale_background_errors(var_scaling3(1), len_scaling3(1), &
 458                                       xb % ds, be3_rf_lengthscale, be % v3)
 459    call da_rescale_background_errors(var_scaling4(1), len_scaling4(1), &
 460                                       xb % ds, be4_rf_lengthscale, be % v4)
 461    call da_rescale_background_errors(var_scaling5(1), len_scaling5(1), &
 462                                       xb % ds, be5_rf_lengthscale, be % v5)
 463
 464    ! 8.0 deallocate input model state:
 465    ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 466
 467    deallocate (be1_eval_loc)
 468    deallocate (be2_eval_loc)
 469    deallocate (be3_eval_loc)
 470    deallocate (be4_eval_loc)
 471    deallocate (be5_eval_loc)
 472
 473    if (vert_corr == vert_corr_2) then
 474       deallocate (be1_eval_glo)
 475       deallocate (be2_eval_glo)
 476       deallocate (be3_eval_glo)
 477       deallocate (be4_eval_glo)
 478       deallocate (be5_eval_glo)
 479
 480       deallocate (be1_evec_loc)
 481       deallocate (be2_evec_loc)
 482       deallocate (be3_evec_loc)
 483       deallocate (be4_evec_loc)
 484       deallocate (be5_evec_loc)
 485
 486       deallocate (be1_evec_glo)
 487       deallocate (be2_evec_glo)
 488       deallocate (be3_evec_glo)
 489       deallocate (be4_evec_glo)
 490       deallocate (be5_evec_glo)
 491
 492    end if
 493
 494    deallocate (bin)
 495    deallocate (bin2d)
 496
 497    if (be % ne > 0) then
 498       be % alpha % mz = be % ne
 499       be % alpha % name = 'alpha'
 500
 501       allocate (be5_eval_loc (1:1,1:num_bins2d))
 502       allocate (be5_eval_glo(1:1))
 503       allocate (be5_evec_loc(1:1,1:1,1:num_bins2d))
 504       allocate (be5_evec_glo(1:1,1:1))
 505
 506       ! Transfer the alpha standard deviation to be % alpha:
 507       allocate (be % alpha % val(1:jy,1:be % ne))
 508       allocate (be % alpha % evec(1:jy,1:kz,1:be % ne))
 509       allocate (be % alpha % rf_alpha(1:be % ne))
 510
 511       be % alpha % val(1:jy,1:be % ne) = alpha_std_dev
 512       be % alpha % evec(1:jy,1:kz,1:be % ne) = 1.0
 513       be % alpha % rf_alpha(1:be % ne) = 1.0
 514
 515       ! Include alpha lengthscale info:
 516       allocate (alpha_rf_lengthscale(1:be % ne))
 517       alpha_rf_lengthscale(:) = 1000.0 * alpha_corr_scale ! Convert km to m.
 518
 519       call da_rescale_background_errors(1.0, 1.0, &
 520          xb % ds, alpha_rf_lengthscale, be % alpha)
 521    else
 522       be % alpha % mz = 0
 523    end if
 524
 525    if (trace_use) call da_trace_exit("da_setup_be_nmm_regional")
 526
 527 end subroutine da_setup_be_nmm_regional
 528
 529