var/gen_be/gen_be_cov2d.f90

   1 program gen_be_cov2d
   2
   3    use da_control, only : stdout, stderr, filename_len
   4    use da_tools_serial, only : da_get_unit,da_advance_cymdh
   5    use da_gen_be, only : da_create_bins
   6
   7    implicit none
   8
   9    character*10        :: start_date, end_date       ! Starting and ending dates.
  10    character*10        :: date, new_date             ! Current date (ccyymmddhh).
  11    character*10        :: variable1                                      ! Variable name
  12    character*10        :: variable2                                      ! Variable name
  13    character(len=filename_len)        :: filename                   ! Input filename.
  14    character*3         :: ce                         ! Member index -> character.
  15    integer             :: ni, nj, nk, nkdum          ! Grid dimensions.
  16    integer             :: i, j, member               ! Loop counters.
  17    integer             :: b                          ! Bin marker.
  18    integer             :: sdate, cdate, edate        ! Starting, current ending dates.
  19    integer             :: interval                   ! Period between dates (hours).
  20    integer             :: ne                         ! Number of ensemble members.
  21    integer             :: bin_type                   ! Type of bin to average over.
  22    integer             :: num_bins                   ! Number of bins (3D fields).
  23    integer             :: num_bins2d                 ! Number of bins (3D fields).
  24    real                :: lat_min, lat_max           ! Used if bin_type = 2 (degrees).
  25    real                :: binwidth_lat               ! Used if bin_type = 2 (degrees).
  26    real                :: hgt_min, hgt_max           ! Used if bin_type = 2 (m).
  27    real                :: binwidth_hgt               ! Used if bin_type = 2 (m).
  28    real                :: coeffa, coeffb             ! Accumulating mean coefficients.
  29    logical             :: first_time                 ! True if first file.
  30
  31    real, allocatable   :: latitude(:,:)              ! Latitude (degrees, from south).
  32    real, allocatable   :: height(:,:,:)              ! Height field.
  33    real, allocatable   :: field1(:,:)                ! Field 1.
  34    real, allocatable   :: field2(:,:)                ! Field 2.
  35    integer, allocatable:: bin(:,:,:)                 ! Bin assigned to each 3D point.
  36    integer, allocatable:: bin2d(:,:)                 ! Bin assigned to each 2D point.
  37    integer, allocatable:: bin_pts2d(:)               ! Number of points in bin (2D fields).
  38    real, allocatable   :: covar(:)                   ! Covariance between input fields.
  39    real, allocatable   :: var(:)                     ! Autocovariance of field.
  40
  41    namelist / gen_be_cov2d_nl / start_date, end_date, interval, &
  42                                 ne, bin_type, &
  43                                 lat_min, lat_max, binwidth_lat, &
  44                                 hgt_min, hgt_max, binwidth_hgt, &
  45                                 variable1, variable2
  46
  47    integer :: ounit,iunit,namelist_unit
  48
  49    stderr = 0
  50    stdout = 6
  51
  52
  53 !---------------------------------------------------------------------------------------------
  54    write(6,'(a)')' [1] Initialize namelist variables and other scalars.'
  55 !---------------------------------------------------------------------------------------------
  56
  57
  58    call da_get_unit(ounit)
  59    call da_get_unit(iunit)
  60    call da_get_unit(namelist_unit)
  61
  62
  63    start_date = '2004030312'
  64    end_date = '2004033112'
  65    interval = 24
  66    ne = 1
  67    bin_type = 5
  68    lat_min = -90.0
  69    lat_max = 90.0
  70    binwidth_lat = 10.0
  71    hgt_min = 0.0
  72    hgt_max = 20000.0
  73    binwidth_hgt = 1000.0
  74    variable1 = 'ps_u'
  75    variable2 = 'ps'
  76
  77    open(unit=namelist_unit, file='gen_be_cov2d_nl.nl', &
  78         form='formatted', status='old', action='read')
  79    read(namelist_unit, gen_be_cov2d_nl)
  80    close(namelist_unit)
  81
  82    read(start_date(1:10), fmt='(i10)')sdate
  83    read(end_date(1:10), fmt='(i10)')edate
  84    write(6,'(4a)')' Computing covariance for fields ', variable1 , ' and ', variable2
  85    write(6,'(4a)') ' Time period is ', start_date, ' to ', end_date
  86    write(6,'(a,i8,a)')' Interval between dates = ', interval, 'hours.'
  87    write(6,'(a,i8)')' Number of ensemble members at each time = ', ne
  88
  89    date = start_date
  90    cdate = sdate
  91
  92 !---------------------------------------------------------------------------------------------
  93    write(6,'(2a)')' [2] Read fields, and calculate correlations'
  94 !---------------------------------------------------------------------------------------------
  95
  96    first_time = .true.
  97
  98    do while ( cdate <= edate )
  99       write(6,'(a,a)')'    Processing data for date ', date
 100
 101       do member = 1, ne
 102
 103          write(ce,'(i3.3)')member
 104
 105 !        Read Full-fields:
 106          filename = 'fullflds/'//date(1:10)
 107          filename = trim(filename)//'.fullflds.e'//ce
 108          open (iunit, file = filename, form='unformatted')
 109
 110          read(iunit)ni, nj, nk
 111          if ( first_time ) then
 112             write(6,'(a,3i8)')'    i, j, k dimensions are ', ni, nj, nk
 113             allocate( latitude(1:ni,1:nj) )
 114             allocate( height(1:ni,1:nj,1:nk) )
 115             allocate( bin(1:ni,1:nj,1:nk) )
 116             allocate( bin2d(1:ni,1:nj) )
 117             allocate( field1(1:ni,1:nj) )
 118             allocate( field2(1:ni,1:nj) )
 119          end if
 120          read(iunit)latitude
 121          read(iunit)height
 122
 123 !        Create and sort into bins:
 124          call da_create_bins( ni, nj, nk, bin_type, num_bins, num_bins2d, bin, bin2d, &
 125                               lat_min, lat_max, binwidth_lat, &
 126                               hgt_min, hgt_max, binwidth_hgt, latitude, height )
 127
 128          close(iunit)
 129
 130          if ( first_time ) then
 131             allocate( bin_pts2d(1:num_bins2d) )
 132             allocate( covar(1:num_bins2d) )
 133             allocate( var(1:num_bins2d) )
 134             bin_pts2d(:) = 0
 135             covar(:) = 0.0
 136             var(:) = 0.0
 137             first_time = .false.
 138          end if
 139
 140 !        Read 2D first field:
 141          filename = trim(variable1)//'/'//date(1:10)
 142          filename = trim(filename)//'.'//trim(variable1)//'.e'//ce//'.01'
 143          open (iunit, file = filename, form='unformatted')
 144          read(iunit)ni, nj, nkdum
 145          read(iunit)field1
 146          close(iunit)
 147
 148 !        Read 2D second field:
 149          filename = trim(variable2)//'/'//date(1:10)
 150          filename = trim(filename)//'.'//trim(variable2)//'.e'//ce//'.01'
 151          open (iunit, file = filename, form='unformatted')
 152          read(iunit)ni, nj, nkdum
 153          read(iunit)field2
 154          close(iunit)
 155
 156 !        Calculate covariances:
 157
 158          do j = 1, nj
 159             do i = 1, ni
 160                b = bin2d(i,j)
 161                bin_pts2d(b) = bin_pts2d(b) + 1
 162                coeffa = 1.0 / real(bin_pts2d(b))
 163                coeffb = real(bin_pts2d(b)-1) * coeffa
 164                covar(b) = coeffb * covar(b) + coeffa * field1(i,j) * field2(i,j)
 165                var(b) = coeffb * var(b) + coeffa * field2(i,j) * field2(i,j)
 166             end do
 167          end do
 168       end do  ! End loop over ensemble members.
 169
 170 !     Calculate next date:
 171       call da_advance_cymdh( date, interval, new_date )
 172       date = new_date
 173       read(date(1:10), fmt='(i10)')cdate
 174    end do     ! End loop over times.
 175
 176    filename = trim(variable1)//'.'//trim(variable2)//'.dat'
 177    open (ounit, file = filename, status='unknown')
 178
 179    do j = 1, nj
 180       do i = 1, ni
 181          b = bin2d(i,j)
 182          write(ounit,'(f16.8)')covar(b) / var(b)
 183       end do
 184    end do
 185
 186 end program gen_be_cov2d