wrfv2_fire/var/gen_be/gen_be_ensmean.f90

   1 program gen_be_ensmean
   2 !
   3 !----------------------------------------------------------------------
   4 !  Purpose: Calculate ensemble mean file from input WRF NETCDF input
   5 !  ensemble members.
   6 !
   7 !  Owner: Dale Barker (NCAR/MMM) - WRF wrappper. Thanks to Cindy Bruyere.
   8 !  Please acknowledge author/institute in work that uses this code.
   9 !
  10 !----------------------------------------------------------------------
  11
  12 #ifdef crayx1
  13 #define iargc ipxfargc
  14 #endif
  15
  16    use da_control, only : stdout, stderr,filename_len
  17    use da_reporting, only : da_error,message
  18
  19    implicit none
  20
  21 #include <netcdf.inc>
  22
  23    integer, parameter    :: max_num_dims = 4          ! Maximum number of dimensions.
  24    integer, parameter    :: max_num_vars = 50         ! Maximum number of variables.
  25    integer, parameter    :: unit = 100                ! Unit number.
  26
  27    character (len=filename_len)   :: directory                 ! General filename stub.
  28    character (len=filename_len)   :: filename                  ! General filename stub.
  29    character (len=filename_len)   :: input_file                ! Input file.
  30    character (len=10)    :: var                       ! Variable to search for.
  31    character (len=3)     :: ce                        ! Member index -> character.
  32
  33    integer               :: num_members               ! Ensemble size.
  34    integer               :: nv                        ! Number of variables.
  35    integer               :: v, member, i              ! Loop counters.
  36    integer               :: length                    ! Filename length.
  37    integer               :: rcode                     ! NETCDF return code.
  38    integer               :: cdfid                     ! NETCDF file IDs.
  39    integer               :: cdfid_mean                ! NETCDF file IDs.
  40    integer               :: cdfid_vari                ! NETCDF file IDs.
  41    integer               :: id_var                    ! NETCDF variable ID.
  42    integer               :: ivtype                    ! 4=integer, 5=real, 6=d.p.
  43    integer               :: ndims                     ! Number of field dimensions.
  44    integer               :: natts                     ! Number of field attributes.
  45    real                  :: member_inv                ! 1 / ensemble size.
  46
  47    character(len=10)     :: cv(1:max_num_vars)        ! Default array of variable names.
  48    integer               :: dimids(1:10)              ! Array of dimension IDs.
  49    integer               :: dims(1:max_num_dims)      ! Array of dimensions.
  50    integer               :: istart(1:max_num_dims)    ! Array of dimension starts.
  51    integer               :: iend(1:max_num_dims)      ! Array of dimension ends.
  52
  53    real (kind=4), allocatable     :: data_r(:,:,:)             ! Data array.
  54    real (kind=4), allocatable     :: data_r_mean(:,:,:)        ! Data array mean.
  55    real (kind=4), allocatable     :: data_r_vari(:,:,:)        ! Data array variance.
  56
  57    namelist / gen_be_ensmean_nl / directory, filename, num_members, nv, cv
  58
  59    stderr = 0
  60    stdout = 6
  61
  62 !---------------------------------------------------------------------------------------------
  63    write(6,'(/a)')' [1] Initialize information.'
  64 !---------------------------------------------------------------------------------------------
  65
  66    directory = './'
  67    filename = 'test'
  68    num_members = 56
  69    nv = 1
  70    cv = "U"
  71
  72    open(unit=unit, file='gen_be_ensmean_nl.nl', &
  73         form='formatted', status='old', action='read')
  74    read(unit, gen_be_ensmean_nl)
  75    close(unit)
  76
  77    write(6,'(a,a)')'   Directory = ', trim(directory)
  78    write(6,'(a,a)')'   filename = ', trim(filename)
  79    write(6,'(a,i4)')'   Number of ensemble members = ', num_members
  80    write(6,'(a,i4)')'   Number of variables to average = ', nv
  81    write(6,'(50a)')'   List of variables to average = ', cv(1:nv)
  82
  83 !  Open template ensemble mean with write access:
  84    input_file = trim(directory)//'/'//trim(filename)
  85    length = len_trim(input_file)
  86    rcode = nf_open(input_file(1:length), NF_WRITE, cdfid_mean )
  87    if ( rcode /= 0) then
  88       write(UNIT=message(1),FMT='(A,A)') &
  89          ' Error opening netcdf file ', input_file(1:length)
  90       call da_error(__FILE__,__LINE__,message(1:1))
  91    end if
  92
  93 !  Open template ensemble variance with write access:
  94    input_file = trim(directory)//'/'//trim(filename)//'.vari'
  95    length = len_trim(input_file)
  96    rcode = nf_open(input_file(1:length), NF_WRITE, cdfid_vari )
  97    if ( rcode /= 0) then
  98       write(UNIT=message(1),FMT='(A,A)') &
  99          ' Error opening netcdf file ', input_file(1:length)
 100       call da_error(__FILE__,__LINE__,message(1:1))
 101    end if
 102
 103 !---------------------------------------------------------------------------------------------
 104    write(6,'(/a)')' [2] Extract necessary fields from WRF ensemble forecasts.'
 105 !---------------------------------------------------------------------------------------------
 106
 107    do v = 1, nv ! Loop over variables to average:
 108       var = cv(v)
 109       write(6,'(2a)')' Computing ensemble mean for variable ', var
 110
 111       do member = 1, num_members
 112          write(UNIT=ce,FMT='(i3.3)')member
 113
 114 !        Open file:
 115          input_file = trim(directory)//'.e'//trim(ce)//'/'//trim(filename)
 116          length = len_trim(input_file)
 117          rcode = nf_open( input_file(1:length), NF_NOWRITE, cdfid )
 118
 119          if ( member == 1 ) then
 120 !           Get variable ID:
 121             rcode = nf_inq_varid ( cdfid, var, id_var )
 122
 123 !           Check variable is in file:
 124             if ( rcode /= 0 ) then
 125                write(UNIT=message(1),FMT='(A,A)') &
 126                   var, ' variable is not in input file'
 127                call da_error(__FILE__,__LINE__,message(1:1))
 128              end if
 129
 130 !            Get dimension information for this field:
 131              dimids = 0
 132              dims = 0
 133              rcode = nf_inq_var( cdfid, id_var, var, ivtype, ndims, dimids, natts )
 134              do i = 1, ndims
 135                 rcode = nf_inq_dimlen( cdfid, dimids(i), dims(i) )
 136              end do
 137              istart = 1
 138              iend = 1
 139              do i = 1, ndims-1
 140                 iend(i) = dims(i)
 141              end do
 142
 143 !            Allocate and initialize data:
 144              if ( ivtype == 5 ) then
 145                 allocate( data_r(iend(1),iend(2),iend(3)))
 146                 allocate( data_r_mean(iend(1),iend(2),iend(3)))
 147                 allocate( data_r_vari(iend(1),iend(2),iend(3)))
 148                 data_r_mean = 0.0
 149                 data_r_vari = 0.0
 150              else
 151                 write(UNIT=message(1),FMT='(A,A)') &
 152                    var, ' variable is not real type'
 153                 call da_error(__FILE__,__LINE__,message(1:1))
 154              end if
 155          end if
 156
 157 !        Calculate accumulating mean and variance:
 158          member_inv = 1.0 / real(member)
 159          call ncvgt( cdfid, id_var, istart, iend, data_r, rcode)
 160          data_r_mean = ( real( member-1 ) * data_r_mean + data_r ) * member_inv
 161          data_r_vari = ( real( member-1 ) * data_r_vari + data_r * data_r ) * member_inv
 162
 163          rcode = nf_close( cdfid )
 164
 165       end do ! member
 166
 167       call ncvpt( cdfid_mean, id_var, istart, iend, data_r_mean, rcode)
 168       call ncvpt( cdfid_vari, id_var, istart, iend, data_r_vari, rcode)
 169       deallocate( data_r )
 170       deallocate( data_r_mean )
 171       deallocate( data_r_vari )
 172
 173    end do ! variable
 174
 175    rcode = nf_close( cdfid_mean )
 176    rcode = nf_close( cdfid_vari )
 177
 178 end program gen_be_ensmean
 179