var/gen_be/gen_be_addmean.f90

   1 program gen_be_addmean
   2 !
   3 !----------------------------------------------------------------------
   4 !  Purpose : add updated ensemble mean and perturbations to get analysis
   5 !
   6 !  Arthur P. Mizzi (NCAR/MMM) April 2011
   7 !
   8 !----------------------------------------------------------------------
   9 ! xlf90 -C -I${NETCDF}/include -L${NETCDF}/lib -lnetcdf -lnetcdff -o da_run_add_mean.exe da_run_add_mean.f90
  10 !
  11    implicit none
  12
  13 #include "netcdf.inc"
  14
  15    integer, parameter    :: max_num_vars = 50              ! Maximum number of variables.
  16    integer, parameter    :: max_num_dims = 20              ! Maximum number of dimensions.
  17    integer, parameter    :: unit = 100                     ! Unit number.
  18
  19    integer               :: num_members                    ! Ensemble size.
  20    integer               :: nv                             ! Number of variables.
  21    integer               :: length                         ! Filename length.
  22    integer               :: rcode                          ! NETCDF return code.
  23    integer               :: cdfid_mean, cdfid_pert         ! NETCDF file IDs.
  24    integer               :: nijkv, nijkv_mean, nijkv_pert  ! Array counter.
  25    integer               :: index                          ! Array index.
  26    integer               :: member, v, o, i, j, k, ijkv    ! Loop counters.
  27    integer               :: ivtype                         ! 4=integer, 5=real, 6=d.p.
  28    integer               :: natts                          ! Number of field attributes.
  29    integer               :: istart_mean(1:max_num_vars)    ! Start index.
  30    integer               :: istart_pert(1:max_num_vars)    ! Start index.
  31    integer               :: id_var_mean(1:max_num_vars)    ! NETCDF variable ID.
  32    integer               :: id_var_pert(1:max_num_vars)    ! NETCDF variable ID.
  33    integer               :: ndims_mean(1:max_num_vars)     ! Number of field dimensions.
  34    integer               :: ndims_pert(1:max_num_vars)     ! Number of field dimensions.
  35    integer               :: one(1:max_num_dims)            ! Array of dimension starts.
  36    integer               :: dimids_mean(1:max_num_dims)    ! Array of dimension IDs.
  37    integer               :: dimids_pert(1:max_num_dims)    ! Array of dimension IDs.
  38    integer               :: dims_mean(1:max_num_vars,1:max_num_dims)      ! Array of dimensions.
  39    integer               :: dims_pert(1:max_num_vars,1:max_num_dims)      ! Array of dimensions.
  40
  41    real (kind=4), allocatable     :: data_r_3d(:,:,:)               ! 3-D Data array.
  42    real (kind=4), allocatable     :: data_r_4d(:,:,:,:)             ! 4-D Data array.
  43    real, allocatable              :: xf_mean(:)                     ! Ensemble mean.
  44    real, allocatable              :: xf_pert(:)                     ! Ensemble perturbations.
  45
  46    character (len=3)     :: ce                        ! Member index -> character.
  47    character (len=10)    :: cv(1:max_num_vars)        ! Default array of variable names.
  48    character (len=160)   :: path                      ! Input file path.
  49    character (len=160)   :: file_mean                 ! Input file_mean name.
  50    character (len=160)   :: file_pert                 ! Input file_pert name.
  51    character (len=160)   :: input_file_mean           ! Input file.
  52    character (len=160)   :: input_file_pert           ! Input file.
  53
  54    namelist / add_mean_nl / num_members, nv, cv, path, file_mean, file_pert
  55 !
  56 ! OPEN AND READ NAMELIST DATA
  57    open(unit=unit, file='add_mean_nl.nl', &
  58         form='formatted', status='old', action='read')
  59    read(unit, add_mean_nl)
  60    close(unit)
  61 !
  62    print *, 'NUM_MEMBERS ', num_members
  63    print *, 'NV ', nv
  64    print *, 'CV ', (trim(cv(i)),i=1,nv)
  65    print *, 'PATH ', trim(path)
  66    print *, 'FILE_MEAN ', trim(file_mean)
  67    print *, 'FILE_PERT ', trim(file_pert)
  68 !
  69 ! OPEN THE ENSEMBLE MEAN FILE
  70    input_file_mean = trim(path)//'/'//trim(file_mean)
  71 !   print *, 'INPUT_FILE_MEAN', trim(input_file_mean)
  72    length = len_trim(input_file_mean)
  73    rcode = nf_open( input_file_mean(1:length), NF_WRITE, cdfid_mean )
  74 !
  75 ! OPEN THE ENSEMBLE PERT FILE
  76    input_file_pert = trim(path)//'/'//trim(file_pert)
  77 !   print *, 'INPUT_FILE_PERT', trim(input_file_pert)
  78    length = len_trim(input_file_pert)
  79    rcode = nf_open( input_file_pert(1:length), NF_WRITE, cdfid_pert )
  80 !
  81 ! LOOP OVER VARIABLES (v) TO GET DIMENSIONS
  82    do v = 1, nv
  83 !
  84 ! GET VARIABLE IDENTIFIER
  85       rcode = nf_inq_varid ( cdfid_mean, cv(v), id_var_mean(v) )
  86       if ( rcode /= 0 ) then
  87          write(UNIT=6,FMT='(A,A)') &
  88          cv(v), ' variable is not in input file mean'
  89          stop
  90       end if
  91       rcode = nf_inq_varid ( cdfid_pert, cv(v), id_var_pert(v) )
  92       if ( rcode /= 0 ) then
  93          write(UNIT=6,FMT='(A,A)') &
  94          cv(v), ' variable is not in input file pert'
  95          stop
  96       end if
  97       if (id_var_mean(v) .ne. id_var_pert(v)) then
  98         print *, 'mean var id not equal pert var id ',v,id_var_mean(v),id_var_pert(v)
  99         stop
 100       end if
 101 !
 102 ! CHECK IF VARIABLE IS TYPE REAL
 103       dimids_mean = 0
 104       rcode = nf_inq_var( cdfid_mean, id_var_mean(v), cv(v), ivtype, ndims_mean(v), dimids_mean, natts )
 105       if ( ivtype /= 5 ) then
 106          write(UNIT=6,FMT='(A,A)') cv(v), ' variable is not real type in file_mean'
 107          stop
 108       end if
 109       dimids_pert = 0
 110       rcode = nf_inq_var( cdfid_pert, id_var_pert(v), cv(v), ivtype, ndims_pert(v), dimids_pert, natts )
 111       if ( ivtype /= 5 ) then
 112          write(UNIT=6,FMT='(A,A)') cv(v), ' variable is not real type in file_pert'
 113          stop
 114       end if
 115 !
 116 ! GET DIMENSIONS FOR THIS VARIABLE
 117       dims_mean(v,:) = 0
 118       do i = 1, ndims_mean(v)
 119          rcode = nf_inq_dimlen( cdfid_mean, dimids_mean(i), dims_mean(v,i) )
 120       end do
 121       dims_pert(v,:) = 0
 122       do i = 1, ndims_pert(v)
 123          rcode = nf_inq_dimlen( cdfid_pert, dimids_pert(i), dims_pert(v,i) )
 124       end do
 125    end do
 126 !
 127 ! SET COLUMN VECTOR INDEXING ARRAY
 128    one = 1
 129    istart_mean(1) = 1
 130    do v = 2, nv
 131       istart_mean(v) = istart_mean(v-1) + product(dims_mean(v-1,1:ndims_mean(v-1)-1))
 132    end do
 133    one = 1
 134    istart_pert(1) = 1
 135    do v = 2, nv
 136       istart_pert(v) = istart_pert(v-1) + product(dims_pert(v-1,1:ndims_pert(v-1)-1))
 137    end do
 138 !
 139 ! CALCULATE SIZE OF COLUMN VECTOR AND ALLOCATE ARRAY
 140    nijkv_mean = istart_mean(nv) + product(dims_mean(nv,1:ndims_mean(nv)-1)) - 1
 141    nijkv_pert = istart_pert(nv) + product(dims_pert(nv,1:ndims_pert(nv)-1)) - 1
 142    nijkv=nijkv_mean
 143    if(nijkv_mean .ne. nijkv_pert) then
 144       print *, 'nijkv_mean not equal nijkv_pert ',nijkv_mean, nijkv_pert
 145       stop
 146    end if
 147
 148    allocate(xf_mean(1:nijkv))
 149    allocate(xf_pert(1:nijkv))
 150 !
 151 ! LOOP OVER VARIABLES TO BE UPDATED
 152    do v = 1, nv
 153       print *, 'process variable ',v
 154 !
 155 ! GET MEAN FIELD VARIABLE IDENTIFIER
 156       rcode = nf_inq_varid ( cdfid_mean, cv(v), id_var_mean(v) )
 157       if ( rcode /= 0 ) then
 158          write(UNIT=6,FMT='(A,A)') &
 159          cv(v), ' variable is not in input file mean'
 160          stop
 161       end if
 162 !      print *, 'save mean field'
 163       index = istart_mean(v)
 164       if(ndims_mean(v) == 3) then
 165          allocate( data_r_3d(dims_mean(v,1),dims_mean(v,2),dims_mean(v,3)))
 166          rcode = nf_get_vara_real( cdfid_mean, id_var_mean(v), one, dims_mean(v,:), data_r_3d)
 167          if(rcode.ne.0) then
 168             print *, 'Error read ',trim(input_file_mean),' variable ',cv(v)
 169             print *, 'rcode ',rcode
 170             stop
 171          end if
 172          do k = 1, dims_mean(v,3)
 173             do j = 1, dims_mean(v,2)
 174                do i = 1, dims_mean(v,1)
 175                   xf_mean(index) = data_r_3d(i,j,k)
 176                   index = index + 1
 177                end do
 178             end do
 179          end do
 180          deallocate( data_r_3d )
 181       end if
 182       if(ndims_mean(v) == 4) then
 183          allocate( data_r_4d(dims_mean(v,1),dims_mean(v,2),dims_mean(v,3),dims_mean(v,4)))
 184          rcode = nf_get_vara_real( cdfid_mean, id_var_mean(v), one, dims_mean(v,:), data_r_4d)
 185          if(rcode.ne.0) then
 186             print *, 'Error read ',trim(input_file_mean),' variable ',cv(v)
 187             print *, 'rcode ',rcode
 188             stop
 189          end if
 190          do k = 1, dims_mean(v,3)
 191             do j = 1, dims_mean(v,2)
 192                do i = 1, dims_mean(v,1)
 193                   xf_mean(index) = data_r_4d(i,j,k,1)
 194                   index = index + 1
 195                end do
 196             end do
 197          end do
 198          deallocate( data_r_4d )
 199       end if
 200 !
 201 ! GET PERT FIELD VARIABLE IDENTIFIER
 202       rcode = nf_inq_varid ( cdfid_pert, cv(v), id_var_pert(v) )
 203       if ( rcode /= 0 ) then
 204          write(UNIT=6,FMT='(A,A)') &
 205          cv(v), ' variable is not in input file pert'
 206          stop
 207       end if
 208 !      print *, 'save pert field'
 209       index = istart_pert(v)
 210       if(ndims_pert(v) == 3) then
 211          allocate( data_r_3d(dims_pert(v,1),dims_pert(v,2),dims_pert(v,3)))
 212          rcode = nf_get_vara_real( cdfid_pert, id_var_pert(v), one, dims_pert(v,:), data_r_3d)
 213          if(rcode.ne.0) then
 214             print *, 'Error read ',trim(input_file_pert),' variable ',cv(v)
 215             print *, 'rcode ',rcode
 216             stop
 217          endif
 218          do k = 1, dims_pert(v,3)
 219             do j = 1, dims_pert(v,2)
 220                do i = 1, dims_pert(v,1)
 221                   xf_pert(index) = data_r_3d(i,j,k)
 222                   index = index + 1
 223                end do
 224             end do
 225          end do
 226          deallocate( data_r_3d )
 227       end if
 228       if(ndims_pert(v) == 4) then
 229          allocate( data_r_4d(dims_pert(v,1),dims_pert(v,2),dims_pert(v,3),dims_pert(v,4)))
 230          rcode = nf_get_vara_real( cdfid_pert, id_var_pert(v), one, dims_pert(v,:), data_r_4d)
 231          if(rcode.ne.0) then
 232             print *, 'Error read ',trim(input_file_pert),' variable ',cv(v)
 233             print *, 'rcode ',rcode
 234             stop
 235          end if
 236          do k = 1, dims_pert(v,3)
 237             do j = 1, dims_pert(v,2)
 238                do i = 1, dims_pert(v,1)
 239                   xf_pert(index) = data_r_4d(i,j,k,1)
 240                   index = index + 1
 241                end do
 242             end do
 243          end do
 244          deallocate( data_r_4d )
 245       end if
 246 !
 247 ! ADD THE MEAN AND PERT FIELDS
 248 !      print *, 'add mean and pert'
 249       do ijkv = 1, nijkv
 250          xf_pert(ijkv) = xf_mean(ijkv) + xf_pert(ijkv)
 251       end do
 252 !
 253 ! WRITE UPDATED FIELD TO FILE_PERT
 254 !      print *, 'write total field'
 255       index = istart_pert(v)
 256       if(ndims_pert(v) == 3) then
 257          allocate( data_r_3d(dims_pert(v,1),dims_pert(v,2),dims_pert(v,3)))
 258          do k = 1, dims_pert(v,3)
 259             do j = 1, dims_pert(v,2)
 260                do i = 1, dims_pert(v,1)
 261                   data_r_3d(i,j,k) = xf_pert(index)
 262                   index = index + 1
 263                end do
 264             end do
 265          end do
 266          call ncvpt( cdfid_pert, id_var_pert(v), one, dims_pert(v,:), data_r_3d, rcode)
 267          deallocate( data_r_3d )
 268       end if
 269       if(ndims_pert(v) == 4) then
 270          allocate( data_r_4d(dims_pert(v,1),dims_pert(v,2),dims_pert(v,3),dims_pert(v,4)))
 271          do k = 1, dims_pert(v,3)
 272             do j = 1, dims_pert(v,2)
 273                do i = 1, dims_pert(v,1)
 274                   data_r_4d(i,j,k,1) = xf_pert(index)
 275                   index = index + 1
 276                end do
 277             end do
 278          end do
 279          call ncvpt( cdfid_pert, id_var_pert(v), one, dims_pert(v,:), data_r_4d, rcode)
 280          deallocate( data_r_4d )
 281       endif
 282    enddo
 283    rcode = nf_close( cdfid_mean )
 284    rcode = nf_close( cdfid_pert )
 285    deallocate(xf_mean)
 286    deallocate(xf_pert)
 287 end program gen_be_addmean