wrfv2_fire/var/da/da_setup_structures/da_setup_be_global.inc

   1 subroutine da_setup_be_global (be)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: Define and allocate components of background errors
   5    !---------------------------------------------------------------------------
   6
   7    implicit none
   8
   9    type (be_type), intent(out) :: be                 ! Back. errors structure
  10
  11    integer                     :: nrec
  12    integer                     :: max_wave_in        ! Dimension of input power
  13    integer                     :: i, j, k, n ! Loop counters
  14    ! real, allocatable   :: height(:,:,:)      ! Height field.
  15    integer, allocatable:: bin(:,:,:)      ! Bin assigned to each 3D point
  16    integer, allocatable:: bin2d(:,:)      ! Bin assigned to each 2D point
  17    integer             :: bin_type        ! Type of bin to average over.
  18    integer             :: num_bins        ! Number of bins (3D fields).
  19    integer             :: num_bins2d      ! Number of bins (3D fields).
  20    logical             :: dummy
  21
  22    real                :: binwidth_lat       ! Used if bin_type = 2 (degrees)
  23    real                :: binwidth_hgt       ! Used if bin_type = 2 (m)
  24    real                :: hgt_min, hgt_max   ! Used if bin_type = 2 (m)
  25    real                :: lat_min, lat_max   ! Used if bin_type = 2 (degrees)
  26
  27    character*10         :: variable
  28    integer              :: ni, nj, nk, b, be_unit
  29    integer, allocatable :: max_wave(:)
  30    real, allocatable    :: evec_g(:,:)
  31    real, allocatable    :: eval_g(:)
  32    real, allocatable    :: evec_loc(:,:,:)
  33    real, allocatable    :: eval_loc(:,:)
  34    real, allocatable    :: regcoeff1(:)
  35    real, allocatable    :: regcoeff2(:,:)
  36    real, allocatable    :: regcoeff3(:,:,:)
  37    real, allocatable    :: power(:)               ! Temporary power spectrum.
  38    real, allocatable    :: power2d(:,:)           ! Temporary power spectrum.
  39
  40    if (trace_use) call da_trace_entry("da_setup_be_global")
  41
  42    call da_message((/"[3.0] Set up background errors (be) for global WRF-Var"/))
  43
  44    be % max_wave = ide / 2 - 1
  45
  46    !---------------------------------------------------------------------
  47    ! [1] Read in be metadata:
  48    !---------------------------------------------------------------------
  49
  50    call da_get_unit(be_unit)
  51    open (unit=be_unit,file="be.dat",status="old",form="unformatted")
  52
  53    read (be_unit, end= 99, err = 100) ni, nj, nk
  54    read (be_unit, err = 100) bin_type
  55    read (be_unit, err = 100) lat_min, lat_max, binwidth_lat
  56    read (be_unit, err = 100) hgt_min, hgt_max, binwidth_hgt
  57    read (be_unit, err = 100) num_bins, num_bins2d
  58
  59    allocate (bin(1:ni,1:nj,1:nk))
  60    allocate (bin2d(1:ni,1:nj))
  61    read (be_unit, err = 100) bin(1:ni,1:nj,1:nk)
  62    read (be_unit, err = 100) bin2d(1:ni,1:nj)
  63
  64    if (ni /= ide .or. nj /= jde .or. nk /= kde) then
  65       call da_error(__FILE__,__LINE__, &
  66          (/"Cannot generate BE at this resolution"/))
  67    end if
  68
  69    !---------------------------------------------------------------------
  70    ! [2] Read in be regression coefficients:
  71    !---------------------------------------------------------------------
  72
  73    allocate (regcoeff1(1:num_bins))
  74    allocate (regcoeff2(1:nk,1:num_bins2d))
  75    allocate (regcoeff3(1:nk,1:nk,1:num_bins2d))
  76
  77    read (be_unit,end= 99,  err = 100) regcoeff1
  78    read (be_unit,end= 99,  err = 100) regcoeff2
  79    read (be_unit,end= 99,  err = 100) regcoeff3
  80
  81    allocate (be%reg_chi(1:jde,1:kde))
  82    allocate (be%reg_ps (1:jde,1:kde))
  83    allocate (be%reg_t  (1:jde,1:kde,1:kde))
  84
  85    ! Fill regression coeff. array
  86    do k=1,kde
  87       do j =1, jde
  88          b = bin(1,j,k) ! Assumes bins averaged over i direction.
  89          be%reg_chi(j,k) = regcoeff1(b)
  90       end do
  91    end do
  92
  93    do j=1,jde
  94       b = bin2d(1,j)
  95       do k=1,kde
  96          be%reg_ps(j,k) = regcoeff2(k,b)
  97       end do
  98    end do
  99
 100    do j=1,jde
 101       b = bin2d(1,j)
 102       do i=1,kde
 103          do k=1,kde
 104             be%reg_t(j,i,k) = regcoeff3(i,k,b)
 105          end do
 106       end do
 107    end do
 108
 109    deallocate(regcoeff1)
 110    deallocate(regcoeff2)
 111    deallocate(regcoeff3)
 112
 113    !---------------------------------------------------------------------
 114    ! [3] Read in be vertical eigenmodes:
 115    !---------------------------------------------------------------------
 116
 117    do nrec = 1, 4
 118       read (be_unit,end= 99,  err = 100) variable
 119       read (be_unit,end= 99,  err = 100) nk, num_bins2d
 120
 121       allocate (evec_g(1:nk,1:nk))
 122       allocate (eval_g(1:nk))
 123       allocate (evec_loc(1:nk,1:nk,num_bins2d))
 124       allocate (eval_loc(1:nk,num_bins2d))
 125
 126       read (be_unit,end= 99, err = 100) evec_g
 127       read (be_unit,end= 99, err = 100) eval_g
 128       read (be_unit,end= 99, err = 100) evec_loc
 129       read (be_unit,end= 99, err = 100) eval_loc
 130
 131       if (nrec == 1) then
 132          be % v1 % name = variable
 133          call da_get_bins_info(nj, nk, bin2d, evec_g, eval_g, &
 134             evec_loc, eval_loc, max_vert_var1, var_scaling1(1), be%v1)
 135
 136       else if (nrec == 2) then
 137          be % v2 % name = variable
 138          call da_get_bins_info(nj, nk, bin2d, evec_g, eval_g, &
 139             evec_loc, eval_loc, max_vert_var2, var_scaling2(1), be%v2)
 140
 141       else if (nrec == 3) then
 142          be % v3 % name = variable
 143          call da_get_bins_info(nj, nk, bin2d, evec_g, eval_g, &
 144             evec_loc, eval_loc, max_vert_var3, var_scaling3(1), be%v3)
 145
 146       else if (nrec == 4) then
 147          be % v4 % name = variable
 148          call da_get_bins_info(nj, nk, bin2d, evec_g, eval_g, &
 149             evec_loc, eval_loc, max_vert_var4, var_scaling4(1), be%v4)
 150       end if
 151
 152       deallocate (evec_g)
 153       deallocate (eval_g)
 154       deallocate (evec_loc)
 155       deallocate (eval_loc)
 156
 157    end do ! loop nrec
 158
 159    deallocate (bin)
 160    deallocate (bin2d)
 161
 162    !---------------------------------------------------------------------
 163    ! [4] Read in be power spectra:
 164    !---------------------------------------------------------------------
 165
 166    allocate(max_wave(1:nk))
 167
 168    do k = 1, nk
 169       read (be_unit) variable
 170       read (be_unit) max_wave_in, nrec
 171       read (be_unit) dummy ! use to preserve file format
 172       if (k == 1) then
 173           allocate (power(0:max_wave_in))                      ! Temporary.
 174           allocate (power2d(0:max_wave_in,1:nk))               ! Temporary.
 175       end if
 176       read (be_unit) power(0:max_wave_in)
 177       power2d(:,k) = power(:)
 178
 179       ! Truncate power spectra:
 180       call da_truncate_spectra(be % max_wave, max_wave_in, power_truncation, &
 181                                 power, max_wave(k))
 182    end do
 183
 184    be % v1 % max_wave = maxval(max_wave(1:nk))
 185    write (unit=stdout,fmt='(/3x,3a,i6)') &
 186       'Horizontal truncation for ', be % v1 % name, ' = ', be % v1 % max_wave
 187    allocate (be % v1 % power(0:be % v1 % max_wave,1:nk))
 188    be % v1 % power(0:be % v1 % max_wave,1:nk) = power2d(0:be % v1 % max_wave,1:nk)
 189    be % v1 % power(0,1:nk) = len_scaling1(1) * be % v1 % power(0,1:nk)
 190
 191    do k = 1, nk
 192       read (be_unit) variable
 193       read (be_unit) max_wave_in, nrec
 194       read (be_unit) dummy ! use to preserve file format
 195       read (be_unit) power(0:max_wave_in)
 196       power2d(:,k) = power(:)
 197
 198       ! Truncate power spectra:
 199       call da_truncate_spectra (be % max_wave, max_wave_in, power_truncation, &
 200                                 power, max_wave(k))
 201    end do
 202
 203    be % v2 % max_wave = maxval(max_wave(1:nk))
 204    write (unit=stdout,fmt='(3x,3a,i6)') &
 205       'Horizontal truncation for ', be % v2 % name, ' = ', be % v2 % max_wave
 206    allocate (be % v2 % power(0:be % v2 % max_wave,1:nk))
 207    be % v2 % power(0:be % v2 % max_wave,1:nk) = power2d(0:be % v2 % max_wave,1:nk)
 208    be % v2 % power(0,1:nk) = len_scaling2(1) * be % v2 % power(0,1:nk)
 209
 210    do k = 1, nk
 211       read (be_unit) variable
 212       read (be_unit) max_wave_in, nrec
 213       read (be_unit) dummy ! use to preserve file format
 214       read (be_unit) power(0:max_wave_in)
 215       power2d(:,k) = power(:)
 216
 217       ! Truncate power spectra:
 218       call da_truncate_spectra (be % max_wave, max_wave_in, power_truncation, &
 219                                 power, max_wave(k))
 220    end do
 221
 222    be % v3 % max_wave = maxval(max_wave(1:nk))
 223    write(unit=stdout,fmt='(3x,3a,i6)') &
 224       'Horizontal truncation for ', be % v3 % name, ' = ', be % v3 % max_wave
 225    allocate (be % v3 % power(0:be % v3 % max_wave,1:nk))
 226    be % v3 % power(0:be % v3 % max_wave,1:nk) = power2d(0:be % v3 % max_wave,1:nk)
 227    be % v3 % power(0,1:nk) = len_scaling3(1) * be % v3 % power(0,1:nk)
 228
 229    do k = 1, nk
 230       read (be_unit) variable
 231       read (be_unit) max_wave_in, nrec
 232       read (be_unit) dummy ! use to preserve file format
 233       read (be_unit) power(0:max_wave_in)
 234       power2d(:,k) = power(:)
 235
 236       ! Truncate power spectra:
 237       call da_truncate_spectra (be % max_wave, max_wave_in, power_truncation, &
 238                                 power, max_wave(k))
 239    end do
 240
 241    be % v4 % max_wave = maxval(max_wave(1:nk))
 242    write (unit=stdout,fmt='(3x,3a,i6)') &
 243       'Horizontal truncation for ', be % v4 % name, ' = ', be % v4 % max_wave
 244    allocate (be % v4 % power(0:be % v4 % max_wave,1:nk))
 245    be % v4 % power(0:be % v4 % max_wave,1:nk) = power2d(0:be % v4 % max_wave,1:nk)
 246    be % v4 % power(0,1:nk) = len_scaling4(1) * be % v4 % power(0,1:nk)
 247
 248    ! ps_u:
 249    read (be_unit) variable
 250    be % v5 % name = variable
 251    be % v5 % mz = 1
 252    if (max_vert_var5 <=  0.0) be % v5 % mz = 0
 253    read (be_unit) max_wave_in, nrec
 254    read (be_unit) dummy ! use to preserve file format
 255    read (be_unit) power(0:max_wave_in)
 256
 257    ! Truncate power spectra:
 258    call da_truncate_spectra (be % max_wave, max_wave_in, power_truncation, &
 259                              power, be % v5 % max_wave)
 260
 261    write (unit=stdout,fmt='(3x,3a,i6)') &
 262       'Horizontal truncation for ', be % v5 % name, ' = ', be % v5 % max_wave
 263    allocate (be % v5 % power(0:be % v5 % max_wave,1))
 264    be % v5 % power(0:be % v5 % max_wave,1) = power(0:be % v5 % max_wave)
 265    be % v5 % power(0,1) = len_scaling5(1) * be%v5%power(0,1)
 266
 267    deallocate(power)
 268
 269    !--------------------------------------------------------------
 270    ! [5] Perform checks on eigenvectors:
 271    !--------------------------------------------------------------
 272
 273    if (test_statistics) then
 274       call da_check_eof_decomposition(be%v1%val_g(:), be%v1%evec_g(:,:),&
 275                                      be%v1%name)
 276       call da_check_eof_decomposition(be%v2%val_g(:), be%v2%evec_g(:,:),&
 277                                      be%v2%name)
 278       call da_check_eof_decomposition(be%v3%val_g(:), be%v3%evec_g(:,:),&
 279                                      be%v3%name)
 280       call da_check_eof_decomposition(be%v4%val_g(:), be%v4%evec_g(:,:),&
 281                                      be%v4%name)
 282    end if
 283
 284    !--------------------------------------------------------------
 285    ! [6] Set up alpha (flow-dependent) control variable "errors":
 286    !--------------------------------------------------------------
 287
 288    if (be % ne > 0) then
 289       be % alpha % mz = be % ne
 290       be % alpha % name = 'alpha'
 291       be % alpha % max_wave = alpha_truncation
 292       if (print_detail_be) then
 293          write(unit=stdout,fmt='(3x,3a,i6)') &
 294             'Horizontal truncation for ', be % alpha % name, ' = ', &
 295             be % alpha % max_wave
 296       end if
 297       allocate (power(0:be % alpha % max_wave))
 298       call da_calc_power_spectrum(be % alpha % max_wave, power)
 299
 300       allocate (be % alpha % power(0:be % alpha % max_wave, be % ne))
 301       do n = 1, be % ne
 302          be % alpha % power(0:be % alpha % max_wave,n) = power(0:be % alpha % max_wave)
 303       end do
 304    end if
 305
 306    deallocate (max_wave)
 307
 308    write (unit=stdout,fmt='(A)') " "
 309
 310    close (be_unit)
 311    call da_free_unit(be_unit)
 312
 313    if (trace_use) call da_trace_exit("da_setup_be_global")
 314
 315    return
 316
 317 99 write (unit=message(1),fmt='(a, i5)')' Unexpected end on BE-unit = ',be_unit
 318    call da_error(__FILE__,__LINE__,message(1:1))
 319
 320 100 write (unit=message(1),fmt='(a, i5)')' Read error on BE-unit = ',be_unit
 321    call da_error(__FILE__,__LINE__,message(1:1))
 322
 323    ! FIX? daft having these here
 324
 325    deallocate(power)
 326    deallocate(power2d)
 327
 328 end subroutine da_setup_be_global
 329
 330