<HTML> <BODY BGCOLOR=#ccccdd LINK=#0000aa VLINK=#0000ff ALINK=#ff0000 ><BASE TARGET="bottom_target"><PRE>
<A NAME='GEN_BE_EP1'><A href='../../html_code/gen_be/gen_be_ep1.f90.html#GEN_BE_EP1' TARGET='top_target'><IMG SRC="../../gif/bar_yellow.gif" border=0></A>
program gen_be_ep1,12
use da_control
, only : stderr,stdout, filename_len
use da_gen_be, only : da_filter_regcoeffs
use da_tools_serial, only : da_get_unit, da_free_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) :: 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 (2D fields).
integer :: num_passes ! Recursive filter passes.
integer :: count ! Counter.
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).
real :: rf_scale ! Recursive filter scale.
real :: count_inv ! 1 / count.
real, allocatable :: psi(:,:,:) ! psi.
real, allocatable :: chi(:,:,:) ! chi.
real, allocatable :: temp(:,:,:) ! Temperature.
real, allocatable :: rh(:,:,:) ! Relative humidity.
real, allocatable :: ps(:,:) ! Surface pressure.
real, allocatable :: chi_u(:,:,:) ! chi.
real, allocatable :: temp_u(:,:,:) ! Temperature.
real, allocatable :: ps_u(:,:) ! Surface pressure.
real, allocatable :: psi_mnsq(:,:,:) ! psi.
real, allocatable :: chi_mnsq(:,:,:) ! chi.
real, allocatable :: temp_mnsq(:,:,:) ! Temperature.
real, allocatable :: rh_mnsq(:,:,:) ! Relative humidity.
real, allocatable :: ps_mnsq(:,:) ! Surface pressure.
real, allocatable :: chi_u_mnsq(:,:,:) ! chi.
real, allocatable :: temp_u_mnsq(:,:,:) ! Temperature.
real, allocatable :: ps_u_mnsq(:,:) ! Surface pressure.
integer, allocatable:: bin(:,:,:) ! Bin assigned to each 3D point.
integer, allocatable:: bin2d(:,:) ! Bin assigned to each 2D point.
real, allocatable :: regcoeff1(:) ! psi/chi regression cooefficient.
real, allocatable :: regcoeff2(:,:) ! psi/ps regression cooefficient.
real, allocatable :: regcoeff3(:,:,:) ! psi/T regression cooefficient.
namelist / gen_be_stage2a_nl / start_date, end_date, interval, &
ne, num_passes, rf_scale
integer :: ounit, gen_be_ounit, namelist_unit, iunit
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(gen_be_ounit)
call da_get_unit(namelist_unit)
start_date = '2004030312'
end_date = '2004033112'
interval = 24
ne = 1
num_passes = 0
rf_scale = 1.0
open(unit=namelist_unit, file='gen_be_stage2a_nl.nl', &
form='formatted', status='old', action='read')
read(namelist_unit, gen_be_stage2a_nl)
close(namelist_unit)
read(start_date(1:10), fmt='(i10)')sdate
read(end_date(1:10), fmt='(i10)')edate
write(6,'(4a)')' Computing control variable fields'
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
date = start_date
cdate = sdate
!---------------------------------------------------------------------------------------------
write(6,'(2a)')' [2] Read regression coefficients and bin information:'
!---------------------------------------------------------------------------------------------
filename = 'be.dat'
open (iunit, file = filename, form='unformatted')
read(iunit)ni, nj, nk
read(iunit)bin_type
read(iunit)lat_min, lat_max, binwidth_lat
read(iunit)hgt_min, hgt_max, binwidth_hgt
read(iunit)num_bins, num_bins2d
allocate( bin(1:ni,1:nj,1:nk) )
allocate( bin2d(1:ni,1:nj) )
allocate( regcoeff1(1:num_bins) )
allocate( regcoeff2(1:nk,1:num_bins2d) )
allocate( regcoeff3(1:nk,1:nk,1:num_bins2d) )
read(iunit)bin(1:ni,1:nj,1:nk)
read(iunit)bin2d(1:ni,1:nj)
read(iunit)regcoeff1
read(iunit)regcoeff2
read(iunit)regcoeff3
close(iunit)
allocate( psi(1:ni,1:nj,1:nk) )
allocate( chi(1:ni,1:nj,1:nk) )
allocate( temp(1:ni,1:nj,1:nk) )
allocate( rh(1:ni,1:nj,1:nk) )
allocate( ps(1:ni,1:nj) )
allocate( chi_u(1:ni,1:nj,1:nk) )
allocate( temp_u(1:ni,1:nj,1:nk) )
allocate( ps_u(1:ni,1:nj) )
if ( num_passes > 0 ) then
write(6,'(a,i4,a)')' [3] Apply ', num_passes, ' pass recursive filter to regression coefficients:'
call da_filter_regcoeffs( ni, nj, nk, num_bins, num_bins2d, num_passes, rf_scale, bin, &
regcoeff1, regcoeff2, regcoeff3 )
else
write(6,'(a)')' [3] num_passes = 0. Bypassing recursive filtering.'
end if
!---------------------------------------------------------------------------------------------
write(6,'(a)')' [4] Read standard fields, and compute unbalanced control variable fields:'
!---------------------------------------------------------------------------------------------
date = start_date
cdate = sdate
do while ( cdate <= edate )
write(6,'(a,a)')' Calculating unbalanced fields for date ', date
do member = 1, ne
write(ce,'(i3.3)')member
! Read psi predictor:
variable = 'psi'
filename = trim(variable)//'/'//date(1:10)
filename = trim(filename)//'.'//trim(variable)//'.e'//ce
open (iunit, file = filename, form='unformatted')
read(iunit)ni, nj, nk
read(iunit)psi
close(iunit)
! Calculate unbalanced chi:
variable = 'chi'
filename = trim(variable)//'/'//date(1:10)
filename = trim(filename)//'.'//trim(variable)//'.e'//ce
open (iunit, file = filename, form='unformatted')
read(iunit)ni, nj, nk
read(iunit)chi
close(iunit)
do k = 1, nk
do j = 1, nj
do i = 1, ni
b = bin(i,j,k)
chi_u(i,j,k) = chi(i,j,k) - regcoeff1(b) * psi(i,j,k)
end do
end do
end do
variable = 'chi_u'
filename = trim(variable)//'/'//date(1:10)
filename = trim(filename)//'.'//trim(variable)//'.e'//ce
open (ounit, file = filename, form='unformatted')
write(ounit)ni, nj, nk
write(ounit)chi_u
close(ounit)
! Calculate unbalanced T:
variable = 't'
filename = trim(variable)//'/'//date(1:10)
filename = trim(filename)//'.'//trim(variable)//'.e'//ce
open (iunit, file = filename, form='unformatted')
read(iunit)ni, nj, nk
read(iunit)temp
close(iunit)
do j = 1, nj
do i = 1, ni
b = bin2d(i,j)
do k = 1, nk
temp_u(i,j,k) = temp(i,j,k) - SUM(regcoeff3(k,1:nk,b) * psi(i,j,1:nk))
end do
end do
end do
variable = 't_u'
filename = trim(variable)//'/'//date(1:10)
filename = trim(filename)//'.'//trim(variable)//'.e'//ce
open (ounit, file = filename, form='unformatted')
write(ounit)ni, nj, nk
write(ounit)temp_u
close(ounit)
! Calculate unbalanced ps:
variable = 'ps'
filename = trim(variable)//'/'//date(1:10)
filename = trim(filename)//'.'//trim(variable)//'.e'//ce//'.01'
open (iunit, file = filename, form='unformatted')
read(iunit)ni, nj, nkdum
read(iunit)ps
close(iunit)
do j = 1, nj
do i = 1, ni
b = bin2d(i,j)
ps_u(i,j) = ps(i,j) - SUM(regcoeff2(1:nk,b) * psi(i,j,1:nk))
end do
end do
variable = 'ps_u'
filename = trim(variable)//'/'//date(1:10)
filename = trim(filename)//'.'//trim(variable)//'.e'//ce//'.01'
open (ounit, file = filename, form='unformatted')
write(ounit)ni, nj, 1
write(ounit)ps_u
close(ounit)
end do ! End loop over ensemble members.
!--------------------------------------------------------------------------------------
! Calculate mean control variables (diagnostics):
!--------------------------------------------------------------------------------------
! Read psi predictor:
variable = 'psi'
filename = trim(variable)//'/'//date(1:10)
filename = trim(filename)//'.'//trim(variable)//'.mean'
open (iunit, file = filename, form='unformatted')
read(iunit)ni, nj, nk
read(iunit)psi
close(iunit)
! Calculate unbalanced chi:
variable = 'chi'
filename = trim(variable)//'/'//date(1:10)
filename = trim(filename)//'.'//trim(variable)//'.mean'
open (iunit, file = filename, form='unformatted')
read(iunit)ni, nj, nk
read(iunit)chi
close(iunit)
do k = 1, nk
do j = 1, nj
do i = 1, ni
b = bin(i,j,k)
chi_u(i,j,k) = chi(i,j,k) - regcoeff1(b) * psi(i,j,k)
end do
end do
end do
variable = 'chi_u'
filename = trim(variable)//'/'//date(1:10)
filename = trim(filename)//'.'//trim(variable)//'.mean'
open (ounit, file = filename, form='unformatted')
write(ounit)ni, nj, nk
write(ounit)chi_u
close(ounit)
! Calculate unbalanced T:
variable = 't'
filename = trim(variable)//'/'//date(1:10)
filename = trim(filename)//'.'//trim(variable)//'.mean'
open (iunit, file = filename, form='unformatted')
read(iunit)ni, nj, nk
read(iunit)temp
close(iunit)
do j = 1, nj
do i = 1, ni
b = bin2d(i,j)
do k = 1, nk
temp_u(i,j,k) = temp(i,j,k) - SUM(regcoeff3(k,1:nk,b) * psi(i,j,1:nk))
end do
end do
end do
variable = 't_u'
filename = trim(variable)//'/'//date(1:10)
filename = trim(filename)//'.'//trim(variable)//'.mean'
open (ounit, file = filename, form='unformatted')
write(ounit)ni, nj, nk
write(ounit)temp_u
close(ounit)
! Calculate unbalanced ps:
variable = 'ps'
filename = trim(variable)//'/'//date(1:10)
filename = trim(filename)//'.'//trim(variable)//'.mean'
open (iunit, file = filename, form='unformatted')
read(iunit)ni, nj, nkdum
read(iunit)ps
close(iunit)
do j = 1, nj
do i = 1, ni
b = bin2d(i,j)
ps_u(i,j) = ps(i,j) - SUM(regcoeff2(1:nk,b) * psi(i,j,1:nk))
end do
end do
variable = 'ps_u'
filename = trim(variable)//'/'//date(1:10)
filename = trim(filename)//'.'//trim(variable)//'.mean'
open (ounit, file = filename, form='unformatted')
write(ounit)ni, nj, 1
write(ounit)ps_u
close(ounit)
! 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.
deallocate( bin )
deallocate( bin2d )
deallocate( regcoeff1 )
deallocate( regcoeff2 )
deallocate( regcoeff3 )
!---------------------------------------------------------------------------------------------
write(6,'(a)')' [5] Compute mean square statistics:'
!---------------------------------------------------------------------------------------------
allocate( psi_mnsq(1:ni,1:nj,1:nk) )
allocate( chi_mnsq(1:ni,1:nj,1:nk) )
allocate( temp_mnsq(1:ni,1:nj,1:nk) )
allocate( rh_mnsq(1:ni,1:nj,1:nk) )
allocate( ps_mnsq(1:ni,1:nj) )
allocate( chi_u_mnsq(1:ni,1:nj,1:nk) )
allocate( temp_u_mnsq(1:ni,1:nj,1:nk) )
allocate( ps_u_mnsq(1:ni,1:nj) )
date = start_date
cdate = sdate
do while ( cdate <= edate )
count = 0
psi_mnsq = 0.0
chi_mnsq = 0.0
temp_mnsq = 0.0
rh_mnsq = 0.0
ps_mnsq = 0.0
chi_u_mnsq = 0.0
temp_u_mnsq = 0.0
ps_u_mnsq = 0.0
do member = 1, ne
write(ce,'(i3.3)')member
count = count + 1
count_inv = 1.0 / real(count)
variable = 'psi'
filename = trim(variable)//'/'//date(1:10)
filename = trim(filename)//'.'//trim(variable)//'.e'//ce
open (iunit, file = filename, form='unformatted')
read(iunit)ni, nj, nk
read(iunit)psi
close(iunit)
variable = 'chi'
filename = trim(variable)//'/'//date(1:10)
filename = trim(filename)//'.'//trim(variable)//'.e'//ce
open (iunit, file = filename, form='unformatted')
read(iunit)ni, nj, nk
read(iunit)chi
close(iunit)
variable = 't'
filename = trim(variable)//'/'//date(1:10)
filename = trim(filename)//'.'//trim(variable)//'.e'//ce
open (iunit, file = filename, form='unformatted')
read(iunit)ni, nj, nk
read(iunit)temp
close(iunit)
variable = 'rh'
filename = trim(variable)//'/'//date(1:10)
filename = trim(filename)//'.'//trim(variable)//'.e'//ce
open (iunit, file = filename, form='unformatted')
read(iunit)ni, nj, nk
read(iunit)rh
close(iunit)
variable = 'ps'
filename = trim(variable)//'/'//date(1:10)
filename = trim(filename)//'.'//trim(variable)//'.e'//ce//'.01'
open (iunit, file = filename, form='unformatted')
read(iunit)ni, nj, nkdum
read(iunit)ps
close(iunit)
variable = 'chi_u'
filename = trim(variable)//'/'//date(1:10)
filename = trim(filename)//'.'//trim(variable)//'.e'//ce
open (iunit, file = filename, form='unformatted')
read(iunit)ni, nj, nk
read(iunit)chi_u
close(iunit)
variable = 't_u'
filename = trim(variable)//'/'//date(1:10)
filename = trim(filename)//'.'//trim(variable)//'.e'//ce
open (iunit, file = filename, form='unformatted')
read(iunit)ni, nj, nk
read(iunit)temp_u
close(iunit)
variable = 'ps_u'
filename = trim(variable)//'/'//date(1:10)
filename = trim(filename)//'.'//trim(variable)//'.e'//ce//'.01'
open (iunit, file = filename, form='unformatted')
read(iunit)ni, nj, nkdum
read(iunit)ps
close(iunit)
! Calculate accumulating mean:
psi_mnsq = ( real( count-1 ) * psi_mnsq + psi * psi ) * count_inv
chi_mnsq = ( real( count-1 ) * chi_mnsq + chi * chi ) * count_inv
temp_mnsq = ( real( count-1 ) * temp_mnsq + temp * temp ) * count_inv
rh_mnsq = ( real( count-1 ) * rh_mnsq + rh * rh ) * count_inv
ps_mnsq = ( real( count-1 ) * ps_mnsq + ps * ps ) * count_inv
chi_u_mnsq = ( real( count-1 ) * chi_u_mnsq + chi_u * chi_u ) * count_inv
temp_u_mnsq = ( real( count-1 ) * temp_u_mnsq + temp_u * temp_u ) * count_inv
ps_u_mnsq = ( real( count-1 ) * ps_u_mnsq + ps_u * ps_u ) * count_inv
end do ! End loop over ensemble members.
psi_mnsq = sqrt(psi_mnsq) ! Convert mnsq to stdv (mean=0)
chi_mnsq = sqrt(chi_mnsq) ! Convert mnsq to stdv (mean=0)
temp_mnsq = sqrt(temp_mnsq) ! Convert mnsq to stdv (mean=0)
rh_mnsq = sqrt(rh_mnsq) ! Convert mnsq to stdv (mean=0)
ps_mnsq = sqrt(ps_mnsq) ! Convert mnsq to stdv (mean=0)
chi_u_mnsq = sqrt(chi_u_mnsq) ! Convert mnsq to stdv (mean=0)
temp_u_mnsq = sqrt(temp_u_mnsq) ! Convert mnsq to stdv (mean=0)
ps_u_mnsq = sqrt(ps_u_mnsq) ! Convert mnsq to stdv (mean=0)
! Write mean square statistics:
filename = 'psi/'//date(1:10)//'.psi.stdv'
open (gen_be_ounit, file = filename, form='unformatted')
write(gen_be_ounit)ni, nj, nk
write(gen_be_ounit)psi_mnsq
close(gen_be_ounit)
filename = 'chi/'//date(1:10)//'.chi.stdv'
open (gen_be_ounit, file = filename, form='unformatted')
write(gen_be_ounit)ni, nj, nk
write(gen_be_ounit)chi_mnsq
close(gen_be_ounit)
filename = 't/'//date(1:10)//'.t.stdv'
open (gen_be_ounit, file = filename, form='unformatted')
write(gen_be_ounit)ni, nj, nk
write(gen_be_ounit)temp_mnsq
close(gen_be_ounit)
filename = 'rh/'//date(1:10)//'.rh.stdv'
open (gen_be_ounit, file = filename, form='unformatted')
write(gen_be_ounit)ni, nj, nk
write(gen_be_ounit)rh_mnsq
close(gen_be_ounit)
filename = 'ps/'//date(1:10)//'.ps.stdv'
open (gen_be_ounit, file = filename, form='unformatted')
write(gen_be_ounit)ni, nj, nk
write(gen_be_ounit)ps_mnsq
close(gen_be_ounit)
filename = 'chi_u/'//date(1:10)//'.chi_u.stdv'
open (gen_be_ounit, file = filename, form='unformatted')
write(gen_be_ounit)ni, nj, nk
write(gen_be_ounit)chi_u_mnsq
close(gen_be_ounit)
filename = 't_u/'//date(1:10)//'.t_u.stdv'
open (gen_be_ounit, file = filename, form='unformatted')
write(gen_be_ounit)ni, nj, nk
write(gen_be_ounit)temp_u_mnsq
close(gen_be_ounit)
filename = 'ps_u/'//date(1:10)//'.ps_u.stdv'
open (gen_be_ounit, file = filename, form='unformatted')
write(gen_be_ounit)ni, nj, nk
write(gen_be_ounit)ps_u_mnsq
close(gen_be_ounit)
! 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.
deallocate( psi )
deallocate( chi )
deallocate( temp )
deallocate( rh )
deallocate( ps )
deallocate( chi_u )
deallocate( temp_u )
deallocate( ps_u )
deallocate( psi_mnsq )
deallocate( chi_mnsq )
deallocate( temp_mnsq )
deallocate( rh_mnsq )
deallocate( ps_mnsq )
deallocate( chi_u_mnsq )
deallocate( temp_u_mnsq )
deallocate( ps_u_mnsq )
call da_free_unit(ounit)
call da_free_unit(iunit)
call da_free_unit(gen_be_ounit)
call da_free_unit(namelist_unit)
end program gen_be_ep1