<HTML> <BODY BGCOLOR=#ccccdd LINK=#0000aa VLINK=#0000ff ALINK=#ff0000 ><BASE TARGET="bottom_target"><PRE>
<A NAME='GEN_BE_HIST'><A href='../../html_code/gen_be/gen_be_hist.f90.html#GEN_BE_HIST' TARGET='top_target'><IMG SRC="../../gif/bar_yellow.gif" border=0></A>
program gen_be_hist,6
!------------------------------------------------------------------------
! Purpose: Computes frequency distributiopn of BE statistics
! Auothor: Syed RH Rizvi (MMM/NESL/NCAR) Date: 02/01/2010
! & Monika Krysta, (CTBTO, Vienna, Austria)
!
! Note: Please acknowledge author/institute in work that uses this code.
!------------------------------------------------------------------------
use da_control
, only : stderr, stdout, filename_len
use da_tools_serial, only : da_get_unit,da_advance_cymdh
implicit none
character*10 :: start_date, end_date ! Starting and ending dates.
character*10 :: date, new_date ! Current date (ccyymmddhh).
character*10 :: variable ! Variable name
character(len=filename_len) :: dat_dir ! Input data directory.
character(len=filename_len) :: filename ! Input filename.
character*3 :: ce ! Member index -> character.
integer :: ni, nj, nk, nkdum ! Grid dimensions.
integer :: i, j, k, member ! Loop counters.
integer :: b ! Bin marker.
integer :: sdate, cdate, edate ! Starting, current ending dates.
integer :: interval ! Period between dates (hours).
integer :: ne ! Number of ensemble members.
integer :: bin_type ! Type of bin to average over.
integer :: num_bins ! Number of bins (3D fields).
integer :: num_bins2d ! Number of bins (3D fields).
real :: lat_min, lat_max ! Used if bin_type = 2 (degrees).
real :: binwidth_lat ! Used if bin_type = 2 (degrees).
real :: hgt_min, hgt_max ! Used if bin_type = 2 (m).
real :: binwidth_hgt ! Used if bin_type = 2 (m).
logical :: first_time ! True if first file.
real, allocatable :: field(:,:,:) ! Field
integer, allocatable:: bin(:,:,:) ! Bin assigned to each 3D point.
integer, allocatable:: bin2d(:,:) ! Bin assigned to each 2D point.
integer, allocatable:: bin_pts(:) ! Number of points in bin (3D fields).
integer, allocatable:: hist(:,:) ! Binned error values
real, allocatable :: var_field(:) ! variance
real, allocatable :: stdev_field(:) ! standard deviation
real, allocatable :: class_hist(:,:) ! Binned error values
integer, allocatable :: rain_class(:,:) ! 2D rain class of perturbation
character*10 :: rainclvar ! Variable name
integer :: Nstdev, N_dim_hist ! Histogram parameters
integer :: intcl ! Histogram loop counter
namelist / gen_be_hist_nl / start_date, end_date, interval, &
ne, variable, Nstdev, N_dim_hist
integer :: ounit,iunit,namelist_unit
stderr = 0
stdout = 6
!---------------------------------------------------------------------------------------------
write(6,'(a)')' [1] Initialize namelist variables and other scalars.'
!---------------------------------------------------------------------------------------------
call da_get_unit(ounit)
call da_get_unit(iunit)
call da_get_unit(namelist_unit)
start_date = '2004030312'
end_date = '2004033112'
interval = 24
ne = 1
variable = 'psi'
dat_dir = '/mmmtmp1/dmbarker'
Nstdev = 5
N_dim_hist = 20
open(unit=namelist_unit, file='gen_be_hist_nl.nl', &
form='formatted', status='old', action='read')
read(namelist_unit, gen_be_hist_nl)
close(namelist_unit)
read(start_date(1:10), fmt='(i10)')sdate
read(end_date(1:10), fmt='(i10)')edate
write(6,'(2a)')' Computing error histogram for field ', variable
write(6,'(4a)') ' Time period is ', start_date, ' to ', end_date
write(6,'(a,i8,a)')' Interval between dates = ', interval, 'hours.'
write(6,'(a,i8)')' Number of ensemble members at each time = ', ne
write(6,'(2(a,i8))')' Parameters of the histogram Nstdev = ', Nstdev,&
' and N_dim_hist = ',N_dim_hist
date = start_date
cdate = sdate
!---------------------------------------------------------------------------------------------
write(6,'(a)')' [2] Accumulate values of errors per bin and vert. level'
!---------------------------------------------------------------------------------------------
first_time = .true.
do while ( cdate <= edate )
write(6,'(a,a)')' Processing data for date ', date
do member = 1, ne
write(ce,'(i3.3)')member
! Read field:
filename = trim(variable)//'/'//date(1:10)
filename = trim(filename)//'.'//trim(variable)//'.e'//ce
open (iunit, file = filename, form='unformatted')
read(iunit)ni, nj, nk
if ( first_time ) then
write(6,'(a,3i8)')' i, j, k dimensions are ', ni, nj, nk
allocate( bin(1:ni,1:nj,1:nk) )
allocate( bin2d(1:ni,1:nj) )
allocate( field(1:ni,1:nj,1:nk) )
! Read bin info:
filename = 'bin.data'
open (iunit+1, file = filename, form='unformatted')
read (iunit+1) bin_type
read (iunit+1) lat_min, lat_max, binwidth_lat
read (iunit+1) hgt_min, hgt_max, binwidth_hgt
read (iunit+1) num_bins, num_bins2d
read (iunit+1) bin(1:ni,1:nj,1:nk)
read (iunit+1) bin2d(1:ni,1:nj)
close(iunit+1)
end if
read(iunit)field
close(iunit)
if ( first_time ) then
write(6,'(a)')' Setup of histogram parameters'
allocate( class_hist(1:nk,1:N_dim_hist) )
allocate( var_field(1:nk) )
allocate( stdev_field(1:nk) )
class_hist = 0.0
var_field = 0.0
stdev_field = 0.0
! finds out global stdev of first field to set up dimensions
do k=1, nk
var_field(k) = 1.0/real(ni*nj-1.0)*sum(field(1:ni,1:nj,k)**2) &
-real(ni*nj)/real(ni*nj-1.0)*(sum(field(1:ni,1:nj,k))/real(ni*nj))**2
stdev_field(k) = sqrt(var_field(k))
write(6,'(a,i4,3(a,e13.5))')' var_field(',k,')= ', var_field(k),&
' lower=',-Nstdev*stdev_field(k),' upper=',Nstdev*stdev_field(k)
do i=1,N_dim_hist
class_hist(k,i)=-Nstdev*stdev_field(k)+2*Nstdev*stdev_field(k)*real(i-1)/real(N_dim_hist-1)
! write(6,'(2(a,i4),a,e13.5') ' class_hist(',k,',',i,')=',class_hist(k,i)
end do
end do
! allocate histogram
allocate( hist(1:N_dim_hist,1:num_bins))
hist(:,:) = 0
!--ym-- 2D rain class of perturbation
if (bin_type==7) allocate( rain_class(1:ni,1:nj) )
!--ym-- 2D rain class of perturbation
first_time = .false.
end if
!--ym-- 2D rain class of perturbation
if (bin_type==7) then
! Read rain_class:
rainclvar = 'raincl'
filename = trim(rainclvar)//'/'//date(1:10)
filename = trim(filename)//'.'//trim(rainclvar)//'.e'//ce//'.01'
open (iunit, file = filename, form='unformatted')
read(iunit)ni, nj, nkdum
read(iunit)rain_class
close(iunit)
! bin depends on rain class of perturbation!
! re-read bin info:
filename = 'bin.data'
open (iunit+1, file = filename, form='unformatted')
read(iunit+1)bin_type
read(iunit+1)lat_min, lat_max, binwidth_lat
read(iunit+1)hgt_min, hgt_max, binwidth_hgt
read(iunit+1)num_bins, num_bins2d
read(iunit+1)bin(1:ni,1:nj,1:nk)
read(iunit+1)bin2d(1:ni,1:nj)
close(iunit+1)
! update it 2D
do j = 1, nj
do i = 1, ni
bin2d(i,j)=rain_class(i,j)*num_bins2d/4+bin2d(i,j)
end do
end do
! update it 3D
do k = 1, nk
do j = 1, nj
do i = 1, ni
bin(i,j,k)=rain_class(i,j)*num_bins/4+bin(i,j,k)
end do
end do
end do
end if
!--ym-- 2D rain class of perturbation
! accumulate values in histogram
do k = 1, nk
do j = 1, nj
do i = 1, ni
b = bin(i,j,k)
! finds out in which bin calls we fall
intcl=int(1+0.5*(N_dim_hist-1.0)*(1+field(i,j,k)/(Nstdev*stdev_field(k))))
if (intcl.ge.1 .and. intcl.le.N_dim_hist) then
hist(intcl,b)=hist(intcl,b)+1
else
! finds out gross errors (to find bug in rh fields)
if (abs(field(i,j,k))>50*Nstdev*stdev_field(k)) then
write(6,'(3(a,i3),2(a,e20.5e3))') ' WARNING Gross error -> err_field(',i," ",j," ",k,&
") =",field(i,j,k)," whereas stdev =",stdev_field(k)
end if
end if
end do
end do
end do
end do ! End loop over ensemble members.
! Calculate next date:
call da_advance_cymdh( date, interval, new_date )
date = new_date
read(date(1:10), fmt='(i10)')cdate
end do ! End loop over times.
!---------------------------------------------------------------------------------------------
write(6,'(a)')' [3] Write out computed histogram'
!---------------------------------------------------------------------------------------------
! write out results
filename = 'hist.'//trim(variable)//'.dat'
open (ounit, file = filename, form='unformatted')
write(ounit)nk,num_bins,N_dim_hist
! split until I do find a working NCL script that reads large binairies
do k=1,nk
write(ounit)class_hist(k,:)
end do
do b=1,num_bins
write(ounit)hist(:,b)
end do
close(ounit)
end program gen_be_hist