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