var/gen_be/gen_be_cov3d3d_contrib.f90

   1 program gen_be_cov3d3d_contrib
   2 !------------------------------------------------------------------------
   3 !  Purpose: Computes contribution of 3D Field on the balanced part
   4 !           of 3D Field for WRFDA cv_options=6
   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    use da_control, only : stderr, stdout, filename_len
  13    use da_tools_serial, only : da_get_unit,da_advance_cymdh
  14
  15    implicit none
  16
  17    character*10        :: start_date, end_date       ! Starting and ending dates.
  18    character*10        :: date, new_date             ! Current date (ccyymmddhh).
  19    character*80        :: variable
  20    character*80        :: regcoeff_var
  21    character*10        :: variable1                  ! Variable name
  22    character*10        :: variable2                  ! Variable name
  23
  24    real, allocatable   :: regcoeff_psi_chi(:)        ! psi/chi    regression cooefficient.
  25    real, allocatable   :: regcoeff_psi_t(:,:,:)      ! psi/t      regression cooefficient.
  26    real, allocatable   :: regcoeff_psi_ps(:,:)       ! psi/ps     regression cooefficient.
  27    real, allocatable   :: regcoeff_psi_rh(:,:,:)     ! psi/rh     regression cooefficient.
  28    real, allocatable   :: regcoeff_chi_u_t(:,:,:)    ! chi_u/t    regression coefficient
  29    real, allocatable   :: regcoeff_chi_u_ps(:,:)     ! chi_u/ps   regression coefficient
  30    real, allocatable   :: regcoeff_chi_u_rh(:,:,:)   ! chi_u/rh   regression coefficient
  31    real, allocatable   :: regcoeff_t_u_rh(:,:,:)     ! t_u/rh     regression coefficient
  32    real, allocatable   :: regcoeff_ps_u_rh(:,:)      ! ps_u/rh    regression coefficient
  33
  34    character(len=filename_len)        :: filename                   ! Input filename.
  35    character*3         :: ce                         ! Member index -> character.
  36    integer             :: ni, nj, nk                 ! Grid dimensions.
  37    integer             :: i, j, k, member            ! Loop counters.
  38    integer             :: b                          ! Bin marker.
  39    integer             :: sdate, cdate, edate        ! Starting, current ending dates.
  40    integer             :: interval                   ! Period between dates (hours).
  41    integer             :: ne                         ! Number of ensemble members.
  42    integer             :: bin_type                   ! Type of bin to average over.
  43    integer             :: num_bins                   ! Number of bins (3D fields).
  44    integer             :: num_bins2d                 ! Number of bins (3D fields).
  45    real                :: lat_min, lat_max           ! Used if bin_type = 2 (degrees).
  46    real                :: binwidth_lat               ! Used if bin_type = 2 (degrees).
  47    real                :: hgt_min, hgt_max           ! Used if bin_type = 2 (m).
  48    real                :: binwidth_hgt               ! Used if bin_type = 2 (m).
  49    real                :: coeffa, coeffb             ! Accumulating mean coefficients.
  50    logical             :: first_time                 ! True if first file.
  51
  52    real, allocatable   :: field1(:,:,:)              ! Field 1.
  53    real, allocatable   :: field2(:,:,:)              ! Field 2 (full)
  54    real, allocatable   :: field2_balanced(:,:,:)   ! Field 2 (unbalanced)
  55    real*4, allocatable :: field(:,:)
  56    integer, allocatable:: bin(:,:,:)                 ! Bin assigned to each 3D point.
  57    integer, allocatable:: bin2d(:,:)                 ! Bin assigned to each 2D point.
  58    integer, allocatable:: bin_pts(:)                 ! Number of points in bin (3D fields).
  59    real, allocatable   :: covar(:)                   ! Covariance between input fields.
  60    real, allocatable   :: var(:)                     ! Autocovariance of field.
  61
  62    real, allocatable   :: regcoeff_field1_field2(:,:,:) ! Regression coefficient
  63
  64    namelist / gen_be_cov3d_nl / start_date, end_date, interval, &
  65                                 ne, variable1, variable2
  66
  67    integer :: ounit,iunit,namelist_unit
  68
  69    stderr = 0
  70    stdout = 6
  71
  72 !---------------------------------------------------------------------------------------------
  73    write(6,'(a)')' [1] Initialize namelist variables and other scalars.'
  74 !---------------------------------------------------------------------------------------------
  75
  76    call da_get_unit(ounit)
  77    call da_get_unit(iunit)
  78    call da_get_unit(namelist_unit)
  79
  80
  81    start_date = '2007070100'
  82    end_date = '2007070500'
  83    interval = 24
  84    ne = 1
  85    variable1 = 'psi'
  86    variable2 = 'chi'
  87
  88    open(unit=namelist_unit, file='gen_be_cov_contrib_nl.nl', &
  89         form='formatted', status='old', action='read')
  90    read(namelist_unit, gen_be_cov3d_nl)
  91    close(namelist_unit)
  92
  93    read(start_date(1:10), fmt='(i10)')sdate
  94    read(end_date(1:10), fmt='(i10)')edate
  95    write(6,'(4a)')' Computing covariance for fields ', variable1 , ' and ', variable2
  96    write(6,'(4a)') ' Time period is ', start_date, ' to ', end_date
  97    write(6,'(a,i8,a)')' Interval between dates = ', interval, 'hours.'
  98    write(6,'(a,i8)')' Number of ensemble members at each time = ', ne
  99
 100
 101    date = start_date
 102    cdate = sdate
 103
 104    member = 1
 105    write(ce,'(i3.3)')member
 106
 107 !  Read grid dimensions from first field file
 108    filename = trim(variable1)//'/'//date(1:10)
 109    filename = trim(filename)//'.'//trim(variable1)//'.e'//ce
 110
 111
 112    open (iunit, file = filename, form='unformatted')
 113    read(iunit)ni, nj, nk
 114    close(iunit)
 115
 116
 117 !  Allocate
 118    write(6,'(a,3i8)')'    i, j, k dimensions are ', ni, nj, nk
 119    allocate( bin(1:ni,1:nj,1:nk) )
 120    allocate( bin2d(1:ni,1:nj) )
 121    allocate( field1(1:ni,1:nj,1:nk) )
 122    allocate( field2(1:ni,1:nj,1:nk) )
 123    allocate( field2_balanced(1:ni,1:nj,1:nk) )
 124    allocate( field(1:ni,1:nj) )
 125
 126 !  Read bin info:
 127    filename = 'bin.data'
 128    open (iunit, file = filename, form='unformatted')
 129    read(iunit)bin_type
 130    read(iunit)lat_min, lat_max, binwidth_lat
 131    read(iunit)hgt_min, hgt_max, binwidth_hgt
 132    read(iunit)num_bins, num_bins2d
 133    read(iunit)bin(1:ni,1:nj,1:nk)
 134    read(iunit)bin2d(1:ni,1:nj)
 135    close(iunit)
 136
 137    allocate( bin_pts(1:num_bins) )
 138    allocate( covar(1:num_bins) )
 139    allocate( var(1:num_bins) )
 140
 141    allocate  (regcoeff_psi_chi(1:num_bins))
 142    allocate  (regcoeff_psi_t(1:nk,1:nk,1:num_bins2d))
 143    allocate  (regcoeff_psi_ps(1:nk,1:num_bins2d))
 144    allocate  (regcoeff_psi_rh(1:nk,1:nk,1:num_bins2d))
 145    allocate  (regcoeff_chi_u_t(1:nk,1:nk,1:num_bins2d))
 146    allocate  (regcoeff_chi_u_ps(1:nk,1:num_bins2d))
 147    allocate  (regcoeff_chi_u_rh(1:nk,1:nk,1:num_bins2d))
 148    allocate  (regcoeff_t_u_rh(1:nk,1:nk,1:num_bins2d))
 149    allocate  (regcoeff_ps_u_rh(1:nk,1:num_bins2d))
 150
 151    allocate( regcoeff_field1_field2(1:nk,1:nk,1:num_bins2d) )
 152
 153    bin_pts(:) = 0
 154    covar(:) = 0.0
 155    var(:) = 0.0
 156
 157 !---------------------------------------------------------------------------------------------
 158    write(6,'(2a)')' [2] Read fields, and calculate correlations'
 159 !---------------------------------------------------------------------------------------------
 160
 161    filename = 'gen_be_stage2.dat'
 162    open(iunit, file = filename, form='unformatted',status='old')
 163
 164    do k = 1, 2
 165    read (iunit)
 166    end do
 167
 168    regcoeff_var='regcoeff_'//trim(variable1)//'_'//trim(variable2)
 169
 170    do k = 1 , 9
 171    read (iunit) variable
 172    select case( trim(adjustl(variable)) )
 173
 174    case ('regcoeff_psi_chi')
 175    read (iunit) regcoeff_psi_chi
 176    case ('regcoeff_psi_t')
 177    read (iunit) regcoeff_psi_t
 178    if(trim(adjustl(variable)) == trim(adjustl(regcoeff_var)) )then
 179    regcoeff_field1_field2 = regcoeff_psi_t
 180    exit
 181    end if
 182    case ('regcoeff_psi_ps')
 183    read (iunit) regcoeff_psi_ps
 184    case ('regcoeff_psi_rh')
 185    read (iunit) regcoeff_psi_rh
 186    if(trim(adjustl(variable)) == trim(adjustl(regcoeff_var)) )then
 187    regcoeff_field1_field2 = regcoeff_psi_rh
 188    exit
 189    end if
 190    case ('regcoeff_chi_u_t')
 191    read (iunit) regcoeff_chi_u_t
 192    if(trim(adjustl(variable)) == trim(adjustl(regcoeff_var)) ) then
 193    regcoeff_field1_field2 = regcoeff_chi_u_t
 194    exit
 195    end if
 196    case ('regcoeff_chi_u_ps')
 197    read (iunit) regcoeff_chi_u_ps
 198    case ('regcoeff_chi_u_rh')
 199    read (iunit) regcoeff_chi_u_rh
 200    if(trim(adjustl(variable)) == trim(adjustl(regcoeff_var)) )then
 201    regcoeff_field1_field2 = regcoeff_chi_u_rh
 202    exit
 203    end if
 204    case ('regcoeff_t_u_rh')
 205    read (iunit) regcoeff_t_u_rh
 206    if(trim(adjustl(variable)) == trim(adjustl(regcoeff_var)) )then
 207    regcoeff_field1_field2 = regcoeff_t_u_rh
 208    exit
 209    end if
 210    case ('regcoeff_ps_u_rh')
 211    read (iunit) regcoeff_ps_u_rh
 212    case default;
 213       write(6,fmt='(a)')'Read problem gen_be_stage2.dat'
 214       write(6,'(a,a)')' Trying to read regression coefficients:',trim(adjustl(variable))
 215       stop
 216    end select
 217    end do
 218
 219    close(iunit)
 220    deallocate (regcoeff_psi_chi)
 221    deallocate (regcoeff_psi_t)
 222    deallocate (regcoeff_psi_ps)
 223    deallocate (regcoeff_psi_rh)
 224    deallocate (regcoeff_chi_u_t)
 225    deallocate (regcoeff_chi_u_ps)
 226    deallocate (regcoeff_chi_u_rh)
 227    deallocate (regcoeff_t_u_rh)
 228    deallocate (regcoeff_ps_u_rh)
 229
 230
 231
 232    do while ( cdate <= edate )
 233
 234       do member = 1, ne
 235
 236          write(ce,'(i3.3)')member
 237
 238 !        Read first field:
 239          filename = trim(variable1)//'/'//date(1:10)
 240          filename = trim(filename)//'.'//trim(variable1)//'.e'//ce
 241          open (iunit, file = filename, form='unformatted')
 242          read(iunit)ni, nj, nk
 243          read(iunit)field1
 244          close(iunit)
 245
 246 !        Read second field:
 247          filename = trim(variable2)//'/'//date(1:10)
 248          filename = trim(filename)//'.'//trim(variable2)//'.e'//ce
 249          open (iunit, file = filename, form='unformatted')
 250          read(iunit)ni, nj, nk
 251          read(iunit)field2
 252          close(iunit)
 253
 254 !        Calculate balanced contribution to field2
 255          do k = 1, nk
 256             do j = 1, nj
 257                do i = 1, ni
 258                   b = bin2d(i,j)
 259                   field2_balanced(i,j,k) = SUM(regcoeff_field1_field2(k,1:nk,b) * field1(i,j,1:nk))
 260                end do
 261             end do
 262          end do
 263
 264 !        Calculate covariances:
 265          do k = 1, nk
 266             do j = 1, nj
 267                do i = 1, ni
 268                   b = bin(i,j,k)
 269                   bin_pts(b) = bin_pts(b) + 1
 270                   coeffa = 1.0 / real(bin_pts(b))
 271                   coeffb = real(bin_pts(b)-1) * coeffa
 272                   covar(b) = coeffb * covar(b) + coeffa * field2_balanced(i,j,k) * field2(i,j,k)
 273                   var(b) = coeffb * var(b) + coeffa * field2(i,j,k) * field2(i,j,k)
 274                end do
 275             end do
 276          end do
 277
 278       end do  ! End loop over ensemble members.
 279
 280 !     Calculate next date:
 281       call da_advance_cymdh( date, interval, new_date )
 282       date = new_date
 283       read(date(1:10), fmt='(i10)')cdate
 284
 285    end do     ! End loop over times.
 286
 287    filename = trim(variable1)//'.'//trim(variable2)//'.bin'
 288    open (ounit, file = filename, form = 'unformatted')
 289    do k = 1, nk
 290       do j = 1, nj
 291          do i = 1, ni
 292             b = bin(i,j,k)
 293             if ( var(b) /= 0.0 ) then
 294                field(i,j) = covar(b) / var(b)
 295             else
 296                field(i,j) = 0.
 297             end if
 298          end do
 299       end do
 300       write(ounit) field
 301    end do
 302    close(ounit)
 303
 304 end program gen_be_cov3d3d_contrib
 305