var/gen_be/gen_be_cov2d3d_contrib.f90

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