var/gen_be/gen_be_diags_read.f90

   1 program gen_be_diags_read
   2 !------------------------------------------------------------------------
   3 !  Purpose: Process background error statistics for display of various
   4 !           diagnostics using WRFDA/graphics/NCL graphics utilities
   5 !
   6 !  Auothor: Syed RH Rizvi (MMM/NESL/NCAR)   Date: 02/01/2010
   7 !           & Monika Krysta, (CTBTO, Vienna, Austria)
   8 !
   9 !  Note: Please acknowledge author/institute in work that uses this code.
  10 !------------------------------------------------------------------------
  11
  12
  13    use da_control, only : do_normalize,filename_len,stderr,stdout,use_rf
  14    use da_tools_serial, only : da_get_unit
  15    use da_gen_be, only : da_print_be_stats_h_regional, &
  16       da_print_be_stats_h_global, da_print_be_stats_p, da_print_be_stats_v
  17
  18    implicit none
  19
  20    character*10        :: variable          ! Variable name
  21    character*80        :: var80c
  22
  23    character*8         :: uh_method         ! Uh_method (power, scale)
  24    integer             :: n_smth_sl                  ! Number of smoothing for scale-length
  25    integer             :: cv_options                 ! WRFDA CV_OPTIONS
  26    character(len=filename_len)        :: filename          ! Input filename.
  27    integer             :: outunit           ! Output unit for diagnostics.
  28    integer             :: ni, nj, nk, nk_3d ! Dimensions read in.
  29    integer             :: bin_type          ! Type of bin to average over. !!!DALE ADD.
  30    integer             :: num_bins          ! Number of 3D bins.
  31    integer             :: num_bins2d        ! Number of 2D bins.
  32    integer             :: k                 ! Loop counter.
  33    integer             :: kdum              ! Dummy vertical index.
  34    integer             :: max_wavenumber    ! Smallest scale required (ni/2 - 1).
  35    real                :: lat_min, lat_max  ! Used if bin_type = 2 (degrees).
  36    real                :: binwidth_lat      ! Used if bin_type = 2 (degrees). !!!DALE ADD..
  37    real                :: binwidth_hgt      ! Used if bin_type = 2 (m). !!!DALE ADD..
  38    real                :: hgt_min, hgt_max  ! Used if bin_type = 2 (m).
  39    real                :: scale_length_ps   ! Scale length for scalar ps or ps_u
  40    logical             :: dummy
  41
  42    integer, allocatable:: bin(:,:,:)        ! Bin assigned to each 3D point.
  43    integer, allocatable:: bin2d(:,:)        ! Bin assigned to each 2D point.
  44    real, allocatable   :: regcoeff_psi_chi(:)        ! psi/chi    regression cooefficient.
  45    real, allocatable   :: regcoeff_psi_t(:,:,:)   ! psi/t   regression cooefficient.
  46    real, allocatable   :: regcoeff_psi_ps(:,:)    ! psi/ps     regression cooefficient.
  47    real, allocatable   :: regcoeff_psi_rh(:,:,:)  ! psi/rh     regression cooefficient.
  48    real, allocatable   :: regcoeff_chi_u_t(:,:,:) ! chi_u/t regression coefficient
  49    real, allocatable   :: regcoeff_chi_u_ps(:,:)  ! chi_u/ps   regression coefficient
  50    real, allocatable   :: regcoeff_chi_u_rh(:,:,:)! chi_u/rh   regression coefficient
  51    real, allocatable   :: regcoeff_t_u_rh(:,:,:)  ! t_u/rh  regression coefficient
  52    real, allocatable   :: regcoeff_ps_u_rh(:,:)   ! ps_u/rh    regression coefficient
  53
  54    real, allocatable   :: e_vec(:,:)        ! Domain-averaged eigenvectors.
  55    real, allocatable   :: e_val(:)          ! Domain-averaged eigenvalues.
  56    real, allocatable   :: e_vec_loc(:,:,:)  ! Latitudinally varying eigenvectors.
  57    real, allocatable   :: e_val_loc(:,:)    ! Latitudinally varying eigenvalues.
  58    real, allocatable   :: total_power(:)    ! Total Power spectrum.
  59    real, allocatable   :: scale_length(:)   ! Scale length for regional application.
  60
  61
  62    namelist / gen_be_diags_nl / cv_options, do_normalize, n_smth_sl, uh_method, use_rf
  63
  64
  65    integer :: iunit,namelist_unit
  66
  67    ! Hardwire because of complicated incrementing of unit numbers writing a file
  68    ! each time
  69    integer, parameter :: ounit = 71
  70
  71    stderr = 0
  72    stdout = 6
  73
  74    uh_method = 'scale'
  75    n_smth_sl = 0
  76    cv_options = 5
  77
  78    call da_get_unit(iunit)
  79    call da_get_unit(namelist_unit)
  80
  81    open(unit=namelist_unit, file='gen_be_diags_nl.nl', &
  82         form='formatted', status='old', action='read')
  83    read(namelist_unit, gen_be_diags_nl)
  84    close(namelist_unit)
  85
  86    filename = 'be.dat'
  87    print '("*** Unit=",i3,3X,"filename=",a40)',iunit, filename
  88    open (iunit, file = filename, form='unformatted')
  89
  90    !----------------------------------------------------------------------------
  91    ! [1] Gather regression coefficients.
  92    !----------------------------------------------------------------------------
  93
  94    ! Read the dimensions:
  95    read(iunit)ni, nj, nk
  96    nk_3d = nk
  97
  98    allocate( bin(1:ni,1:nj,1:nk) )
  99    allocate( bin2d(1:ni,1:nj) )
 100
 101    ! Read bin info:
 102
 103    read(iunit)bin_type
 104    read(iunit)lat_min, lat_max, binwidth_lat
 105    read(iunit)hgt_min, hgt_max, binwidth_hgt
 106    read(iunit)num_bins, num_bins2d
 107    read(iunit)bin(1:ni,1:nj,1:nk)
 108    read(iunit)bin2d(1:ni,1:nj)
 109
 110
 111    if ( cv_options /= 7 ) then
 112       ! 1.1 Read in regression coefficients
 113       allocate  (regcoeff_psi_chi(1:num_bins))
 114       allocate  (regcoeff_psi_t(1:nk,1:nk,1:num_bins2d))
 115       allocate  (regcoeff_psi_ps(1:nk,1:num_bins2d))
 116       allocate  (regcoeff_psi_rh(1:nk,1:nk,1:num_bins2d))
 117       allocate  (regcoeff_chi_u_t(1:nk,1:nk,1:num_bins2d))
 118       allocate  (regcoeff_chi_u_ps(1:nk,1:num_bins2d))
 119       allocate  (regcoeff_chi_u_rh(1:nk,1:nk,1:num_bins2d))
 120       allocate  (regcoeff_t_u_rh(1:nk,1:nk,1:num_bins2d))
 121       allocate  (regcoeff_ps_u_rh(1:nk,1:num_bins2d))
 122
 123       regcoeff_psi_chi = 0.
 124       regcoeff_psi_t   = 0.
 125       regcoeff_psi_ps  = 0.
 126       regcoeff_psi_rh  = 0.
 127       regcoeff_chi_u_t = 0.
 128       regcoeff_chi_u_ps= 0.
 129       regcoeff_chi_u_rh= 0.
 130       regcoeff_t_u_rh  = 0.
 131       regcoeff_ps_u_rh = 0.
 132
 133       if ( cv_options == 5 ) then
 134          read (iunit) regcoeff_psi_chi
 135          read (iunit) regcoeff_psi_ps
 136          read (iunit) regcoeff_psi_t
 137       else
 138          do k = 1 , 9
 139             read (iunit) var80c
 140             select case( trim(adjustl(var80c)) )
 141             case ('regcoeff_psi_chi')
 142                read (iunit) regcoeff_psi_chi
 143             case ('regcoeff_psi_t')
 144                read (iunit) regcoeff_psi_t
 145             case ('regcoeff_psi_ps')
 146                read (iunit) regcoeff_psi_ps
 147             case ('regcoeff_psi_rh')
 148                read (iunit) regcoeff_psi_rh
 149             case ('regcoeff_chi_u_t')
 150                read (iunit) regcoeff_chi_u_t
 151             case ('regcoeff_chi_u_ps')
 152                read (iunit) regcoeff_chi_u_ps
 153             case ('regcoeff_chi_u_rh')
 154                read (iunit) regcoeff_chi_u_rh
 155             case ('regcoeff_t_u_rh')
 156                read (iunit) regcoeff_t_u_rh
 157             case ('regcoeff_ps_u_rh')
 158                read (iunit) regcoeff_ps_u_rh
 159             case default;
 160                write(6,fmt='(A,A)')'Error in reading regression coefficients for: ',trim(adjustl(var80c))
 161                stop
 162             end select
 163          end do
 164       end if
 165
 166       outunit = ounit + 100
 167       call da_print_be_stats_p( outunit, ni, nj, nk, num_bins, num_bins2d, &
 168                     bin, bin2d, regcoeff_psi_chi,regcoeff_psi_ps,regcoeff_psi_t)
 169    else
 170       !outunit is advanced by 3 within da_print_be_stats_p; since we don't call
 171       !that for cv_options=7 we need to account for it
 172       outunit = ounit + 103
 173    end if
 174
 175    !----------------------------------------------------------------------------
 176    ! [2] Gather vertical error eigenvectors, eigenvalues.
 177    !----------------------------------------------------------------------------
 178
 179    read(iunit)variable
 180    read(iunit)nk, num_bins2d
 181
 182    allocate( e_vec(1:nk,1:nk) )
 183    allocate( e_val(1:nk) )
 184    allocate( e_vec_loc(1:nk,1:nk,1:num_bins2d) )
 185    allocate( e_val_loc(1:nk,1:num_bins2d) )
 186
 187    read(iunit)e_vec
 188    read(iunit)e_val
 189    read(iunit)e_vec_loc
 190    read(iunit)e_val_loc
 191    call da_print_be_stats_v( outunit, variable, nk, num_bins2d, &
 192                              e_vec, e_val, e_vec_loc, e_val_loc )
 193    read(iunit)variable
 194    read(iunit)nk, num_bins2d
 195    read(iunit)e_vec
 196    read(iunit)e_val
 197    read(iunit)e_vec_loc
 198    read(iunit)e_val_loc
 199    call da_print_be_stats_v( outunit, variable, nk, num_bins2d, &
 200                              e_vec, e_val, e_vec_loc, e_val_loc )
 201
 202    read(iunit)variable
 203    read(iunit)nk, num_bins2d
 204    read(iunit)e_vec
 205    read(iunit)e_val
 206    read(iunit)e_vec_loc
 207    read(iunit)e_val_loc
 208    call da_print_be_stats_v( outunit, variable, nk, num_bins2d, &
 209                              e_vec, e_val, e_vec_loc, e_val_loc )
 210
 211    read(iunit)variable
 212    read(iunit)nk, num_bins2d
 213    read(iunit)e_vec
 214    read(iunit)e_val
 215    read(iunit)e_vec_loc
 216    read(iunit)e_val_loc
 217    call da_print_be_stats_v( outunit, variable, nk, num_bins2d, &
 218                              e_vec, e_val, e_vec_loc, e_val_loc )
 219
 220    deallocate( e_vec )
 221    deallocate( e_val )
 222    deallocate( e_vec_loc )
 223    deallocate( e_val_loc )
 224
 225    if (uh_method /= 'power') then
 226       ! 2d fields: ps_u, ps:
 227
 228       read(iunit)variable
 229       read(iunit)nk, num_bins2d
 230
 231       allocate( e_vec(1:nk,1:nk) )
 232       allocate( e_val(1:nk) )
 233       allocate( e_vec_loc(1:nk,1:nk,1:num_bins2d) )
 234       allocate( e_val_loc(1:nk,1:num_bins2d) )
 235
 236       read(iunit)e_vec
 237       read(iunit)e_val
 238       read(iunit)e_vec_loc
 239       read(iunit)e_val_loc
 240       call da_print_be_stats_v( outunit, variable, nk, num_bins2d, &
 241                                 e_vec, e_val, e_vec_loc, e_val_loc )
 242       deallocate( e_vec )
 243       deallocate( e_val )
 244       deallocate( e_vec_loc )
 245       deallocate( e_val_loc )
 246    end if
 247    ! To assign the dimension nk for 3d fields:
 248    nk = nk_3d
 249
 250    if (uh_method == 'power') then
 251        write(6,'(/a)') '[3] Gather horizontal error power spectra.'
 252
 253       do k = 1, nk
 254          read(iunit)variable
 255          read(iunit)max_wavenumber, kdum
 256          if ( k == 1 ) allocate( total_power(0:max_wavenumber) )
 257          read(iunit) dummy ! to preserve namelist format
 258          read(iunit)total_power(:)
 259          call da_print_be_stats_h_global( outunit, variable, k, max_wavenumber, total_power )
 260       end do
 261
 262       do k = 1, nk
 263          read(iunit)variable
 264          read(iunit)max_wavenumber, kdum
 265          read(iunit) dummy ! to preserve namelist format
 266          read(iunit)total_power(:)
 267          call da_print_be_stats_h_global( outunit, variable, k, max_wavenumber, total_power )
 268       end do
 269
 270       do k = 1, nk
 271          read(iunit)variable
 272          read(iunit)max_wavenumber, kdum
 273          read(iunit) dummy ! to preserve namelist format
 274          read(iunit)total_power(:)
 275          call da_print_be_stats_h_global( outunit, variable, k, max_wavenumber, total_power )
 276       end do
 277
 278       do k = 1, nk
 279          read(iunit)variable
 280          read(iunit)max_wavenumber, kdum
 281          read(iunit) dummy ! to preserve namelist format
 282          read(iunit)total_power(:)
 283          call da_print_be_stats_h_global( outunit, variable, k, max_wavenumber, total_power )
 284       end do
 285
 286       read(iunit)variable
 287       read(iunit)max_wavenumber, kdum
 288       read(iunit) dummy ! to preserve namelist format
 289       read(iunit)total_power(:)
 290       call da_print_be_stats_h_global( outunit, variable, k, max_wavenumber, total_power )
 291
 292    else if (uh_method == 'scale   ') then
 293
 294       allocate (scale_length(1:nk))
 295
 296       ! psi/u:
 297       read(iunit) variable
 298       read(iunit) scale_length
 299       call da_print_be_stats_h_regional( outunit, variable, nk, scale_length )
 300
 301       ! chi_u/v:
 302       read(iunit) variable
 303       read(iunit) scale_length
 304       call da_print_be_stats_h_regional( outunit, variable, nk, scale_length )
 305
 306       ! t_u/t:
 307       read(iunit) variable
 308       read(iunit) scale_length
 309       call da_print_be_stats_h_regional( outunit, variable, nk, scale_length )
 310
 311       ! rh:
 312       read(iunit) variable
 313       read(iunit) scale_length
 314       call da_print_be_stats_h_regional( outunit, variable, nk, scale_length )
 315
 316       ! ps_u/ps:
 317       read(iunit) variable
 318       read(iunit) scale_length_ps ! This is a scalar, can be ps_u or ps depending on cv_options
 319       write(6,'(3a,i5)')' Scale length for variable ', trim(variable), ' in unit ', outunit
 320       write(outunit,'(a,i4,1pe15.5)')trim(variable), 1, scale_length_ps
 321       outunit = outunit + 1
 322       write(6,*)
 323
 324       deallocate (scale_length)
 325
 326    else
 327       write(6,'(a,": uh_method=",a," not implemented.")')__FILE__,uh_method
 328    endif
 329
 330    close(iunit)
 331
 332 end program gen_be_diags_read