var/convertor/average_be/average_be.F

   1 program average_be
   2
   3    use be_type
   4    use readwrf_module
   5
   6    implicit none
   7
   8    integer :: i, j, k, jj, kk, n, b, istatus
   9    integer :: mapproj, west_east_dim, south_north_dim, idim, jdim, kdim, idim_stag, jdim_stag
  10    integer :: nmax
  11    integer, dimension(1) :: nmax_a
  12    real :: d
  13    real :: avg_psfc
  14    real :: delta_lat
  15    real, allocatable, dimension(:) :: wt
  16    real, allocatable, dimension(:) :: temp_stats_eval_loc, temp_xb_eval_loc
  17    real, allocatable, dimension(:,:) :: temp_stats_evec_loc, temp_xb_evec_loc
  18    real, allocatable, dimension(:) :: stats_regcoeff1, stats_regcoeff2
  19    real, allocatable, dimension(:) :: xb_regcoeff1, xb_regcoeff2
  20    real, allocatable, dimension(:,:) :: stats_regcoeff3, xb_regcoeff3
  21    real*4 :: dx, dy, cen_lat, cen_lon, stand_lon, true1, true2, ratio, miycors, mjxcors
  22    real, allocatable, dimension(:) :: xb_pres, stats_pres
  23    character (len=256) :: be_fname
  24    type (be_dat) :: new_h_be, new_be
  25    type (be_dat), allocatable, dimension(:) :: be
  26
  27    real  :: lat1,coslat1
  28
  29    ! Namelist variables
  30    integer :: nbins
  31    character (len=256) :: fg_file_name
  32
  33    namelist /nl_average_be/ nbins, fg_file_name
  34
  35
  36 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  37 ! Read namelist
  38 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  39    open(10,file='namelist.average_be',status='old')
  40    read(10,nl_average_be)
  41    close(10)
  42
  43
  44 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  45 ! Read in BE stats for each latitude bin
  46 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  47    allocate(be(nbins))
  48    do i=1,nbins
  49       write(6,*) '***** Reading in BE stats for bin ',i
  50       write(be_fname,'(a,i1)') 'be.dat.',i
  51       call rd_be_cv_5(be_fname, be(i))
  52    end do
  53
  54
  55 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  56 ! Read in a first guess file to get information about the model domain
  57 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  58    write(6,*) '***** Reading in first guess file ', trim(fg_file_name)
  59    istatus = readwrf(fg_file_name, west_east_dim, south_north_dim, dx, dy, cen_lat, cen_lon, &
  60                  stand_lon, true1, true2, mapproj, idim, jdim, kdim, idim_stag, jdim_stag, &
  61                  ratio, miycors, mjxcors)
  62
  63
  64 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  65 ! Set up BE structure for interpolated statistics
  66 !
  67 ! new_h_be will hold the horizontally interpolated BE statistics
  68 ! new_be   will hold the vertically and horizontally interpolated BE statistics
  69 !
  70 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  71    write(6,*) '***** Allocating BE structures'
  72    new_be%ni = idim
  73    new_be%nj = jdim
  74    new_be%nk = kdim
  75    new_be%nk_2d = 1
  76
  77    new_be%bin_type = 5
  78    new_be%num_bins = kdim
  79    new_be%num_bins2d = 1
  80
  81    new_be%lat_min = minval(wrf_xlat)   ! Minimum latitude in wrfinput file
  82    new_be%lat_max =  maxval(wrf_xlat)  ! Maximum latitude in wrfinput file
  83    new_be%binwidth_lat = (new_be%lat_max - new_be%lat_min) / real(jdim-1)
  84
  85    new_be%hgt_min = 0.0
  86    new_be%hgt_max =  0.0
  87    new_be%binwidth_hgt = 0.0
  88
  89    new_be%variable(:) = be(1)%variable(:)
  90
  91    allocate( new_be%bin (idim,jdim,kdim) )
  92    allocate( new_be%bin2d (idim,jdim) )
  93
  94    allocate( new_be%regcoeff1 (new_be%num_bins) )
  95    allocate( new_be%regcoeff2 (kdim,new_be%num_bins2d) )
  96    allocate( new_be%regcoeff3 (kdim,kdim,new_be%num_bins2d) )
  97
  98    allocate ( new_be%be1_eval_loc (1:jdim,1:kdim) )
  99    allocate ( new_be%be2_eval_loc (1:jdim,1:kdim) )
 100    allocate ( new_be%be3_eval_loc (1:jdim,1:kdim) )
 101    allocate ( new_be%be4_eval_loc (1:jdim,1:kdim) )
 102    allocate ( new_be%be5_eval_loc (1:jdim,1:1 ) )
 103
 104    allocate ( new_be%be1_eval_glo (1:kdim) )
 105    allocate ( new_be%be2_eval_glo (1:kdim) )
 106    allocate ( new_be%be3_eval_glo (1:kdim) )
 107    allocate ( new_be%be4_eval_glo (1:kdim) )
 108    allocate ( new_be%be5_eval_glo (1:1) )
 109
 110    allocate ( new_be%be1_evec_loc (1:jdim,1:kdim,1:kdim))
 111    allocate ( new_be%be2_evec_loc (1:jdim,1:kdim,1:kdim))
 112    allocate ( new_be%be3_evec_loc (1:jdim,1:kdim,1:kdim))
 113    allocate ( new_be%be4_evec_loc (1:jdim,1:kdim,1:kdim))
 114    allocate ( new_be%be5_evec_loc (1:jdim,1: 1,1: 1))
 115
 116    allocate ( new_be%be1_evec_glo (1:kdim,1:kdim) )
 117    allocate ( new_be%be2_evec_glo (1:kdim,1:kdim) )
 118    allocate ( new_be%be3_evec_glo (1:kdim,1:kdim) )
 119    allocate ( new_be%be4_evec_glo (1:kdim,1:kdim) )
 120    allocate ( new_be%be5_evec_glo (1:1,1:1) )
 121
 122    allocate ( new_be%be1_rf_lengthscale (1:kdim) )
 123    allocate ( new_be%be2_rf_lengthscale (1:kdim) )
 124    allocate ( new_be%be3_rf_lengthscale (1:kdim) )
 125    allocate ( new_be%be4_rf_lengthscale (1:kdim) )
 126    allocate ( new_be%be5_rf_lengthscale (1:kdim) )
 127
 128
 129
 130    allocate( new_h_be%bin (idim,jdim,be(1)%nk) )
 131    allocate( new_h_be%bin2d (idim,jdim) )
 132
 133    ! NB: The regression coefficient arrays of new_h_be are never actually used
 134    allocate( new_h_be%regcoeff1 (new_be%num_bins) )
 135    allocate( new_h_be%regcoeff2 (be(1)%nk,new_be%num_bins2d) )
 136    allocate( new_h_be%regcoeff3 (be(1)%nk,be(1)%nk,new_be%num_bins2d) )
 137
 138    allocate ( new_h_be%be1_eval_loc (1:jdim,1:be(1)%nk) )
 139    allocate ( new_h_be%be2_eval_loc (1:jdim,1:be(1)%nk) )
 140    allocate ( new_h_be%be3_eval_loc (1:jdim,1:be(1)%nk) )
 141    allocate ( new_h_be%be4_eval_loc (1:jdim,1:be(1)%nk) )
 142    allocate ( new_h_be%be5_eval_loc (1:jdim,1:1 ) )
 143
 144    allocate ( new_h_be%be1_eval_glo (1:be(1)%nk) )
 145    allocate ( new_h_be%be2_eval_glo (1:be(1)%nk) )
 146    allocate ( new_h_be%be3_eval_glo (1:be(1)%nk) )
 147    allocate ( new_h_be%be4_eval_glo (1:be(1)%nk) )
 148    allocate ( new_h_be%be5_eval_glo (1:1) )
 149
 150    allocate ( new_h_be%be1_evec_loc (1:jdim,1:be(1)%nk,1:be(1)%nk))
 151    allocate ( new_h_be%be2_evec_loc (1:jdim,1:be(1)%nk,1:be(1)%nk))
 152    allocate ( new_h_be%be3_evec_loc (1:jdim,1:be(1)%nk,1:be(1)%nk))
 153    allocate ( new_h_be%be4_evec_loc (1:jdim,1:be(1)%nk,1:be(1)%nk))
 154    allocate ( new_h_be%be5_evec_loc (1:jdim,1: 1,1: 1))
 155
 156    allocate ( new_h_be%be1_evec_glo (1:be(1)%nk,1:be(1)%nk) )
 157    allocate ( new_h_be%be2_evec_glo (1:be(1)%nk,1:be(1)%nk) )
 158    allocate ( new_h_be%be3_evec_glo (1:be(1)%nk,1:be(1)%nk) )
 159    allocate ( new_h_be%be4_evec_glo (1:be(1)%nk,1:be(1)%nk) )
 160    allocate ( new_h_be%be5_evec_glo (1:1,1:1) )
 161
 162    allocate ( new_h_be%be1_rf_lengthscale (1:be(1)%nk) )
 163    allocate ( new_h_be%be2_rf_lengthscale (1:be(1)%nk) )
 164    allocate ( new_h_be%be3_rf_lengthscale (1:be(1)%nk) )
 165    allocate ( new_h_be%be4_rf_lengthscale (1:be(1)%nk) )
 166    allocate ( new_h_be%be5_rf_lengthscale (1:be(1)%nk) )
 167
 168
 169 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 170 ! Horizontally interpolate BE statistics
 171 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 172    write(6,*) '***** Horizontally interpolating statistics'
 173
 174    !
 175    ! Determine what percent of the regional domain falls within each global bin
 176    !
 177    allocate(wt(nbins))
 178    wt(:) = 0.0
 179    do i=1,idim
 180       do j=1,jdim
 181          do n=1,nbins
 182             if (wrf_xlat(i,j) >= (real(n-1)*be(n)%binwidth_lat + be(n)%lat_min) .and. &
 183                 wrf_xlat(i,j) <= (real(n)*be(n)%binwidth_lat + be(n)%lat_min)) then
 184                wt(n) = wt(n) + 1.0
 185             end if
 186          end do
 187       end do
 188    end do
 189    do n=1,nbins
 190       wt(n) = wt(n) / real(idim*jdim)
 191       write(0,*) wt(n),' percent of regional domain is within band ', n
 192    end do
 193
 194    ! Set nmax to the global bin that contains largest part of regional domain
 195    nmax_a = maxloc(wt)
 196    nmax = nmax_a(1)
 197
 198
 199    ! Loop over original dimension of evals and evecs
 200    do k=1,be(1)%nk
 201
 202       !
 203       ! Interpolate global eigenvectors
 204       !
 205
 206       ! Loop over eigenvector component
 207       do kk=1,be(1)%nk
 208          new_h_be%be1_evec_glo(k,kk) = be(nmax)%be1_evec_glo(k,kk)
 209          new_h_be%be2_evec_glo(k,kk) = be(nmax)%be2_evec_glo(k,kk)
 210          new_h_be%be3_evec_glo(k,kk) = be(nmax)%be3_evec_glo(k,kk)
 211          new_h_be%be4_evec_glo(k,kk) = be(nmax)%be4_evec_glo(k,kk)
 212
 213          if (k == 1 .and. kk == 1) then
 214             new_h_be%be5_evec_glo(1,1) = be(nmax)%be5_evec_glo(1,1)
 215          end if
 216       end do
 217
 218       !
 219       ! Interpolate global eigenvalues
 220       !
 221
 222       new_h_be%be1_eval_glo(k) = be(nmax)%be1_eval_glo(k)
 223       new_h_be%be2_eval_glo(k) = be(nmax)%be2_eval_glo(k)
 224       new_h_be%be3_eval_glo(k) = be(nmax)%be3_eval_glo(k)
 225       new_h_be%be4_eval_glo(k) = be(nmax)%be4_eval_glo(k)
 226       if (k == 1) new_h_be%be5_eval_glo(k) = be(nmax)%be5_eval_glo(k)
 227
 228       !
 229       ! Interpolate lengthscales
 230       !
 231       new_h_be%be1_rf_lengthscale(k) = 0.0
 232       new_h_be%be2_rf_lengthscale(k) = 0.0
 233       new_h_be%be3_rf_lengthscale(k) = 0.0
 234       new_h_be%be4_rf_lengthscale(k) = 0.0
 235
 236       do n=1,nbins
 237          lat1 = (real(n)-0.5)*be(n)%binwidth_lat + be(n)%lat_min
 238          coslat1 = cos(lat1*deg_to_rad)
 239
 240          new_h_be%be1_rf_lengthscale(k) = new_h_be%be1_rf_lengthscale(k) + coslat1*be(n)%be1_rf_lengthscale(k) * wt(n)
 241          new_h_be%be2_rf_lengthscale(k) = new_h_be%be2_rf_lengthscale(k) + coslat1*be(n)%be2_rf_lengthscale(k) * wt(n)
 242          new_h_be%be3_rf_lengthscale(k) = new_h_be%be3_rf_lengthscale(k) + coslat1*be(n)%be3_rf_lengthscale(k) * wt(n)
 243          new_h_be%be4_rf_lengthscale(k) = new_h_be%be4_rf_lengthscale(k) + coslat1*be(n)%be4_rf_lengthscale(k) * wt(n)
 244          if (k == 1) &
 245          new_h_be%be5_rf_lengthscale(k) = new_h_be%be5_rf_lengthscale(k) + coslat1*be(n)%be5_rf_lengthscale(k) * wt(n)
 246       end do
 247    end do
 248
 249    !
 250    ! Interpolate local eigenvectors
 251    !
 252    do j=1,new_be%nj   ! new j latitudes
 253
 254       ! Set d to the row in the global domain that is closest in latitude to current xb j index
 255       d = (wrf_xlat(1,j)-be(1)%lat_min)/(be(1)%lat_max-be(1)%lat_min)*real(be(1)%nj)
 256
 257       do k=1,be(1)%nk       ! original dimension of evals and evecs
 258       do kk=1,be(1)%nk   ! kk is eigenvector component index
 259
 260          new_h_be%be1_evec_loc(j,k,kk) = be(1)%be1_evec_loc(nint(d),k,kk)
 261          new_h_be%be2_evec_loc(j,k,kk) = be(1)%be2_evec_loc(nint(d),k,kk)
 262          new_h_be%be3_evec_loc(j,k,kk) = be(1)%be3_evec_loc(nint(d),k,kk)
 263          new_h_be%be4_evec_loc(j,k,kk) = be(1)%be4_evec_loc(nint(d),k,kk)
 264          if (k == 1 .and. kk == 1) new_h_be%be5_evec_loc(j,k,kk) = be(1)%be5_evec_loc(nint(d),k,kk)
 265
 266       end do
 267       end do
 268    end do
 269
 270    !
 271    ! Interpolate local eigenvalues
 272    !
 273    do j=1,new_be%nj   ! new j latitudes
 274       d = (wrf_xlat(1,j)-be(1)%lat_min)/(be(1)%lat_max-be(1)%lat_min)*real(be(1)%nj)
 275
 276       do k=1,be(1)%nk       ! original dimension of evals and evecs
 277
 278          new_h_be%be1_eval_loc(j,k) = be(1)%be1_eval_loc(nint(d),k)
 279          new_h_be%be2_eval_loc(j,k) = be(1)%be2_eval_loc(nint(d),k)
 280          new_h_be%be3_eval_loc(j,k) = be(1)%be3_eval_loc(nint(d),k)
 281          new_h_be%be4_eval_loc(j,k) = be(1)%be4_eval_loc(nint(d),k)
 282          if (k == 1) new_h_be%be5_eval_loc(j,k) = be(1)%be5_eval_loc(nint(d),k)
 283
 284       end do
 285    end do
 286
 287    !
 288    ! Set up bins
 289    !
 290     do k=1,kdim
 291        new_be%bin(:,:,k) = k
 292     end do
 293     new_be%bin2d(:,:) = 1
 294
 295
 296    write(6,*) '** Interpolating stats for domain dimensioned ',idim,jdim,kdim
 297    write(6,*) 'Domain lower-left at (lat,lon)=',wrf_xlat(1,1),wrf_xlong(1,1)
 298    write(6,*) 'Domain upper-right at (lat,lon)=',wrf_xlat(idim,jdim),wrf_xlong(idim,jdim)
 299    write(6,*) 'znu from model = ', znu
 300    write(6,*) 'model ptop = ', ptop
 301
 302
 303 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 304 ! Vertically interpolate BE statistics
 305 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 306    write(6,*) '***** Vertically interpolating statistics'
 307
 308    !
 309    ! Fill in array with mean pressure on "xb" levels
 310    !
 311    allocate(xb_pres(kdim))
 312    avg_psfc = sum(wrf_psfc(:,:))/real(idim*jdim)
 313    do k=1,kdim
 314       xb_pres(k) = (avg_psfc - ptop)*znu(k) + ptop
 315       write(6,*) ' xb_pres(',k,')= ',xb_pres(k)
 316    end do
 317
 318    !
 319    ! Fill in array with mean pressure on stats levels
 320    !
 321    allocate(stats_pres(be(1)%nk))
 322    do k=1,be(1)%nk
 323       stats_pres(k) = (be(1)%psfc - be(1)%ptop)*be(1)%znu(k) + be(1)%ptop
 324       write(6,*) ' stats_pres(',k,')= ',stats_pres(k)
 325    end do
 326
 327
 328    !
 329    ! Interpolate global eigenvectors, global eigenvalues, and lengthscales
 330    !
 331    call da_interpolate_stats( be(1)%nk,                    &   ! Number of levels in stats.
 332                               kdim,                        &   ! Number of levels in xb.
 333                               stats_pres,                  &   ! Mean pressure on stats levs.
 334                               xb_pres,                     &   ! Mean pressure on xb levs.
 335                               new_h_be%be1_evec_glo,       &   ! Eigenvectors of vert B.
 336                               new_h_be%be1_eval_glo,       &   ! Eigenvalues of vert B.
 337                               new_h_be%be1_rf_lengthscale, &   ! Correlation scale.
 338                               new_be%be1_evec_glo,         &   ! New eigenvectors.
 339                               new_be%be1_eval_glo,         &   ! New eigenvalues.
 340                               new_be%be1_rf_lengthscale,   &   ! New lengthscales.
 341                               dx                           )   ! grid distance in m
 342
 343    call da_interpolate_stats( be(1)%nk,                    &   ! Number of levels in stats.
 344                               kdim,                        &   ! Number of levels in xb.
 345                               stats_pres,                  &   ! Mean pressure on stats levs.
 346                               xb_pres,                     &   ! Mean pressure on xb levs.
 347                               new_h_be%be2_evec_glo,       &   ! Eigenvectors of vert B.
 348                               new_h_be%be2_eval_glo,       &   ! Eigenvalues of vert B.
 349                               new_h_be%be2_rf_lengthscale, &   ! Correlation scale.
 350                               new_be%be2_evec_glo,         &   ! New eigenvectors.
 351                               new_be%be2_eval_glo,         &   ! New eigenvalues.
 352                               new_be%be2_rf_lengthscale,   &   ! New lengthscales.
 353                               dx                           )   ! grid distance in m
 354
 355    call da_interpolate_stats( be(1)%nk,                    &   ! Number of levels in stats.
 356                               kdim,                        &   ! Number of levels in xb.
 357                               stats_pres,                  &   ! Mean pressure on stats levs.
 358                               xb_pres,                     &   ! Mean pressure on xb levs.
 359                               new_h_be%be3_evec_glo,       &   ! Eigenvectors of vert B.
 360                               new_h_be%be3_eval_glo,       &   ! Eigenvalues of vert B.
 361                               new_h_be%be3_rf_lengthscale, &   ! Correlation scale.
 362                               new_be%be3_evec_glo,         &   ! New eigenvectors.
 363                               new_be%be3_eval_glo,         &   ! New eigenvalues.
 364                               new_be%be3_rf_lengthscale,   &   ! New lengthscales.
 365                               dx                           )   ! grid distance in m
 366
 367    call da_interpolate_stats( be(1)%nk,                    &   ! Number of levels in stats.
 368                               kdim,                        &   ! Number of levels in xb.
 369                               stats_pres,                  &   ! Mean pressure on stats levs.
 370                               xb_pres,                     &   ! Mean pressure on xb levs.
 371                               new_h_be%be4_evec_glo,       &   ! Eigenvectors of vert B.
 372                               new_h_be%be4_eval_glo,       &   ! Eigenvalues of vert B.
 373                               new_h_be%be4_rf_lengthscale, &   ! Correlation scale.
 374                               new_be%be4_evec_glo,         &   ! New eigenvectors.
 375                               new_be%be4_eval_glo,         &   ! New eigenvalues.
 376                               new_be%be4_rf_lengthscale,   &   ! New lengthscales.
 377                               dx                           )   ! grid distance in m
 378
 379    new_be%be5_rf_lengthscale(1) = new_h_be%be5_rf_lengthscale(1)
 380    do k=2,kdim
 381       new_be%be5_rf_lengthscale(k) = 0.0
 382    end do
 383    new_be%be5_evec_glo(:,:) = new_h_be%be5_evec_glo(:,:)
 384    new_be%be5_eval_glo(:)   = new_h_be%be5_eval_glo(:)
 385
 386
 387    !
 388    ! Interpolate local eigenvectors, and local eigenvalues
 389    !
 390    allocate(temp_stats_evec_loc(be(1)%nk,be(1)%nk))
 391    allocate(temp_stats_eval_loc(be(1)%nk))
 392    allocate(temp_xb_evec_loc(kdim,kdim))
 393    allocate(temp_xb_eval_loc(kdim))
 394
 395    do j=1,jdim
 396       temp_stats_evec_loc(:,:) = new_h_be%be1_evec_loc(j,:,:)
 397       temp_stats_eval_loc(:)   = new_h_be%be1_eval_loc(j,:)
 398       call da_interpolate_stats( be(1)%nk,                     &   ! Number of levels in stats.
 399                                  kdim,                         &   ! Number of levels in xb.
 400                                  stats_pres,                   &   ! Mean pressure on stats levs.
 401                                  xb_pres,                      &   ! Mean pressure on xb levs.
 402                                  temp_stats_evec_loc,          &   ! Eigenvectors of vert B.
 403                                  temp_stats_eval_loc,          &   ! Eigenvalues of vert B.
 404                                  new_h_be%be1_rf_lengthscale,  &   ! Correlation scale.
 405                                  temp_xb_evec_loc,             &   ! New eigenvectors.
 406                                  temp_xb_eval_loc,             &   ! New eigenvalues.
 407                                  new_be%be1_rf_lengthscale,    &   ! New lengthscales.
 408                                  dx                            )   ! grid distance in m
 409       new_be%be1_evec_loc(j,:,:) = temp_xb_evec_loc(:,:)
 410       new_be%be1_eval_loc(j,:)   = temp_xb_eval_loc(:)
 411
 412       temp_stats_evec_loc(:,:) = new_h_be%be2_evec_loc(j,:,:)
 413       temp_stats_eval_loc(:)   = new_h_be%be2_eval_loc(j,:)
 414       call da_interpolate_stats( be(1)%nk,                     &   ! Number of levels in stats.
 415                                  kdim,                         &   ! Number of levels in xb.
 416                                  stats_pres,                   &   ! Mean pressure on stats levs.
 417                                  xb_pres,                      &   ! Mean pressure on xb levs.
 418                                  temp_stats_evec_loc,          &   ! Eigenvectors of vert B.
 419                                  temp_stats_eval_loc,          &   ! Eigenvalues of vert B.
 420                                  new_h_be%be2_rf_lengthscale,  &   ! Correlation scale.
 421                                  temp_xb_evec_loc,             &   ! New eigenvectors.
 422                                  temp_xb_eval_loc,             &   ! New eigenvalues.
 423                                  new_be%be2_rf_lengthscale,    &   ! New lengthscales.
 424                                  dx                            )   ! grid distance in m
 425       new_be%be2_evec_loc(j,:,:) = temp_xb_evec_loc(:,:)
 426       new_be%be2_eval_loc(j,:)   = temp_xb_eval_loc(:)
 427
 428       temp_stats_evec_loc(:,:) = new_h_be%be3_evec_loc(j,:,:)
 429       temp_stats_eval_loc(:)   = new_h_be%be3_eval_loc(j,:)
 430       call da_interpolate_stats( be(1)%nk,                     &   ! Number of levels in stats.
 431                                  kdim,                         &   ! Number of levels in xb.
 432                                  stats_pres,                   &   ! Mean pressure on stats levs.
 433                                  xb_pres,                      &   ! Mean pressure on xb levs.
 434                                  temp_stats_evec_loc,          &   ! Eigenvectors of vert B.
 435                                  temp_stats_eval_loc,          &   ! Eigenvalues of vert B.
 436                                  new_h_be%be3_rf_lengthscale,  &   ! Correlation scale.
 437                                  temp_xb_evec_loc,             &   ! New eigenvectors.
 438                                  temp_xb_eval_loc,             &   ! New eigenvalues.
 439                                  new_be%be3_rf_lengthscale,    &   ! New lengthscales.
 440                                  dx                            )   ! grid distance in m
 441       new_be%be3_evec_loc(j,:,:) = temp_xb_evec_loc(:,:)
 442       new_be%be3_eval_loc(j,:)   = temp_xb_eval_loc(:)
 443
 444       temp_stats_evec_loc(:,:) = new_h_be%be4_evec_loc(j,:,:)
 445       temp_stats_eval_loc(:)   = new_h_be%be4_eval_loc(j,:)
 446       call da_interpolate_stats( be(1)%nk,                     &   ! Number of levels in stats.
 447                                  kdim,                         &   ! Number of levels in xb.
 448                                  stats_pres,                   &   ! Mean pressure on stats levs.
 449                                  xb_pres,                      &   ! Mean pressure on xb levs.
 450                                  temp_stats_evec_loc,          &   ! Eigenvectors of vert B.
 451                                  temp_stats_eval_loc,          &   ! Eigenvalues of vert B.
 452                                  new_h_be%be4_rf_lengthscale,  &   ! Correlation scale.
 453                                  temp_xb_evec_loc,             &   ! New eigenvectors.
 454                                  temp_xb_eval_loc,             &   ! New eigenvalues.
 455                                  new_be%be4_rf_lengthscale,    &   ! New lengthscales.
 456                                  dx                            )   ! grid distance in m
 457       new_be%be4_evec_loc(j,:,:) = temp_xb_evec_loc(:,:)
 458       new_be%be4_eval_loc(j,:)   = temp_xb_eval_loc(:)
 459
 460    end do
 461
 462    new_be%be5_evec_loc(:,:,:) = new_h_be%be5_evec_loc(:,:,:)
 463    new_be%be5_eval_loc(:,:)   = new_h_be%be5_eval_loc(:,:)
 464
 465    allocate( stats_regcoeff1(1:be(1)%nk) )
 466    allocate( stats_regcoeff2(1:be(1)%nk) )
 467    allocate( stats_regcoeff3(1:be(1)%nk, 1:be(1)%nk) )
 468
 469    stats_regcoeff1 = 0.
 470    stats_regcoeff2 = 0.
 471    stats_regcoeff3 = 0.
 472    do k=1,be(1)%nk
 473       do n=1,nbins
 474          kk = n + (k-1) * nbins
 475          stats_regcoeff1(k)= stats_regcoeff1(k) + be(1)%regcoeff1(kk)*wt(n)
 476          stats_regcoeff2(k)= stats_regcoeff2(k) + be(1)%regcoeff2(k,n)*wt(n)
 477          do i=1,be(1)%nk
 478             stats_regcoeff3(i,k)= stats_regcoeff3(i,k) + be(1)%regcoeff3(i,k,n)*wt(n)
 479          end do
 480       end do
 481    end do
 482
 483    allocate( xb_regcoeff1(1:new_be%nk) )
 484    allocate( xb_regcoeff2(1:new_be%nk) )
 485    allocate( xb_regcoeff3(1:new_be%nk, 1:new_be%nk) )
 486
 487    call da_interpolate_regcoeff( be(1)%nk, new_be%nk, stats_pres, xb_pres,    &
 488                                  stats_regcoeff1, stats_regcoeff2, stats_regcoeff3, &
 489                                  xb_regcoeff1, xb_regcoeff2, xb_regcoeff3 )
 490
 491 ! Now transfer regression coefficients to new_be array
 492
 493    do k=1,new_be%nk
 494       new_be%regcoeff1 (k)  = xb_regcoeff1(k)
 495       new_be%regcoeff2 (k,:)= xb_regcoeff2(k)
 496       do kk=1,new_be%nk
 497          new_be%regcoeff3(k,kk,:)= xb_regcoeff3(k,kk)
 498       end do
 499    end do
 500
 501    deallocate(temp_stats_evec_loc)
 502    deallocate(temp_stats_eval_loc)
 503    deallocate(temp_xb_evec_loc)
 504    deallocate(temp_xb_eval_loc)
 505
 506
 507 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 508 ! Write out interpolated be statistics to new file
 509 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 510    be_fname = 'be.dat'
 511    write(6,*) '***** Writing interpolated BE stats to ', trim(be_fname)
 512    call wr_be_cv_5(be_fname, new_be)
 513
 514
 515 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 516 ! Clean up
 517 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 518    do i=1,nbins
 519       call free_be_dat(be(i))
 520    end do
 521    deallocate(be)
 522    call free_be_dat(new_be)
 523    call free_be_dat(new_h_be)
 524    deallocate(wt)
 525    deallocate(xb_pres)
 526    deallocate(stats_pres)
 527
 528    stop
 529
 530 end program average_be