data_io.f90

Go to the documentation of this file.

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! Time-stamp: <2016-10-18 15:23:32 pbrowne>
!!!
!!!    Collection of subroutines to deal with i/o
!!!    Copyright (C) 2014  Philip A. Browne
!!!
!!!    This program is free software: you can redistribute it and/or modify
!!!    it under the terms of the GNU General Public License as published by
!!!    the Free Software Foundation, either version 3 of the License, or
!!!    (at your option) any later version.
!!!
!!!    This program is distributed in the hope that it will be useful,
!!!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!!!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!!!    GNU General Public License for more details.
!!!
!!!    You should have received a copy of the GNU General Public License
!!!    along with this program.  If not, see <http://www.gnu.org/licenses/>.
!!!
!!!    Email: p.browne @ reading.ac.uk
!!!    Mail:  School of Mathematical and Physical Sciences,
!!!           University of Reading,
!!!       Reading, UK
!!!       RG6 6BB
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine default_get_observation_data(y,t)
  use output_empire, only : unit_obs,emp_e
  use timestep_data
  use pf_control
  use sizes
  implicit none
  integer, parameter :: rk = kind(1.0d0)
  integer, intent(in) :: t
  real(kind=rk), dimension(obs_dim), intent(out) :: y
  integer :: ios
  character(14) :: filename

  write(filename,'(A,i7.7)') 'obs_ts_',tsdata%next_ob_timestep

  open(unit_obs,file=filename,iostat=ios,action='read',status='old',form='unformatted')
  if(ios .ne. 0)  then
     write(emp_e,*) 'PARTICLE FILTER DATA ERROR!!!!! Cannot open file ',filename
     write(emp_e,*) 'Check it exists. I need a lie down.'
     stop
  end if
  read(unit_obs) y
  close(unit_obs)
end subroutine default_get_observation_data

subroutine save_observation_data(y)
  use output_empire, only : unit_obs,emp_e
  use pf_control
  use sizes
  implicit none
  integer, parameter :: rk = kind(1.0d0)
  real(kind=rk), dimension(obs_dim), intent(in) :: y
  integer :: ios
  character(14) :: filename

  write(filename,'(A,i7.7)') 'obs_ts_',pf%timestep

  open(unit_obs,file=filename,iostat=ios,action='write',status='replace',form='unformatted')
  if(ios .ne. 0)  then
     write(emp_e,*) 'PARTICLE FILTER DATA ERROR!!!!! Cannot open file ',filename
     write(emp_e,*) 'Very strange that I couldnt open it. Im going to stop now.'
     stop
  end if
  write(unit_obs) y
  close(unit_obs)

end subroutine save_observation_data

subroutine get_truth(x)
  use output_empire, only : unit_truth,emp_e
  use timestep_data
  use pf_control
  use sizes
  implicit none
  integer, parameter :: rk = kind(1.0d0)
  real(kind=rk), dimension(state_dim), intent(out) :: x
  integer :: ios
  if(pf%timestep .eq. 0) then
     print*,'opening pf_truth'
     open(unit_truth,file='pf_truth',iostat=ios,action='read')
     if(ios .ne. 0)  then
        write(emp_e,*) 'PARTICLE FILTER DATA ERROR!!!!! Cannot open file pf_truth'
        write(emp_e,*) 'Very strange that I couldnt open it. Im going to stop now.'
        stop
     end if
  end if
  read(unit_truth,*) x
  call flush(unit_truth)
  !if(pf%timestep .eq. pf%time_obs*pf%time_bwn_obs) then
  if(tsdata%completed_timesteps .eq. tsdata%total_timesteps) then
     close(unit_truth)
     print*,'closing pf_truth'
  end if
end subroutine get_truth

subroutine save_truth(x)
  use output_empire, only : unit_truth,unit_mean,emp_e
  use timestep_data
  use pf_control
  use sizes
  implicit none
  integer, parameter :: rk = kind(1.0d0)
  real(kind=rk), dimension(state_dim), intent(in) :: x
  integer :: ios
  if(pf%timestep .eq. 0) then
     print*,'opening pf_truth'
     open(unit_truth,file='pf_truth',iostat=ios,action='write',status='replace')
     if(ios .ne. 0)  then
        write(emp_e,*) 'PARTICLE FILTER DATA ERROR!!!!! Cannot open file pf_truth'
        write(emp_e,*) 'Very strange that I couldnt open it. Im going to stop now.'
        stop
     end if
  end if
  write(unit_truth,*) x
  call flush(unit_truth)
  !if(pf%timestep .eq. pf%time_obs*pf%time_bwn_obs) then
  if(tsdata%completed_timesteps .eq. tsdata%total_timesteps) then
     close(unit_truth)
     print*,'closing pf_truth'
  end if
end subroutine save_truth



subroutine output_from_pf
  use output_empire, only : unit_weight,unit_mean,emp_e
  use timestep_data
  use matrix_pf
  use pf_control
  use sizes
  use comms
  implicit none
  include 'mpif.h'
  real(kind=kind(1.0D0)), dimension(state_dim) :: mean,mtemp
  integer :: ios,particle,mpi_err
  character(20) :: filename


  if(pf%timestep .eq. 0) then

     if(pf%output_weights) then
        write(filename,'(A,i2.2)') 'ensemble_weights_',pfrank
        open(unit_weight,file=filename,iostat=ios,action='write',status='replace')
        if(ios .ne. 0)  then
           write(emp_e,*) 'PARTICLE FILTER DATA ERROR!!!!! Cannot open file pf_out'
           write(emp_e,*) 'Very strange that I couldnt open it. Im going to stop now.'
           stop
        end if
     end if !end if output_weights
     call read_matrix_pf_information
  end if

  if(pf%output_weights) then
     write(unit_weight,'(i6.6,A)',advance='no') pf%timestep,' '
     do ios = 1,pf%count-1
        write(unit_weight,'(i6.6,A,e21.15,A)',advance='no') pf%particles(ios),' ',pf&
             &%weight(pf%particles(ios)),' '
     end do
     write(unit_weight,'(i6.6,A,e21.15)',advance='yes') pf%particles(pf%count),' ',pf&
          &%weight(pf%particles(pf%count))
     call flush(unit_weight)

     if(tsdata%completed_timesteps .eq. tsdata%total_timesteps) close(unit_weight)  
  end if !end if output_weights
     
  if(pf%use_mean .and. pf_ens_rank .eq. 0) then
     if(pf%timestep .eq. 0) then
        open(unit_mean,file='pf_mean',iostat=ios,action='write',status='replace')
        if(ios .ne. 0)  then
           write(emp_e,*) 'PARTICLE FILTER DATA ERROR!!!!! Cannot open file pf_mean'
           write(emp_e,*) 'Very strange that I couldnt open it. Im going to stop now.'
           stop
        end if
     end if
  end if

  if(pf%use_mean .or. (pf%use_spatial_rmse .and. .not. pf%gen_data)) then
     mtemp = 0.0d0
     do particle = 1,pf%count
        !mean(:) = mean(:) + pf%psi(:,particle)*exp(-pf&
        !     &%weight(particles(particle)))
        mtemp(:) = mtemp(:) + pf%psi(:,particle)/real(pf%nens)
     end do

     call mpi_allreduce(mtemp,mean,state_dim,mpi_double_precision,mpi_sum&
             &,pf_ens_comm,mpi_err)
    
  end if

  if(pf%use_mean .and. pf_ens_rank .eq. 0) then
     write(unit_mean,*) mean(:)
     call flush(unit_mean)
     if(tsdata%completed_timesteps .eq. tsdata%total_timesteps) close(unit_mean)
  end if

  if(pf_ens_rank .eq. 0 .and. pf%use_spatial_rmse .and. .not. pf%gen_data) call output_spatial_rmse(mean)

  if(pf%use_mean .and. pf%use_variance) call output_variance(mean)

  if(comm_version .ne. 3) then !empire_v3 models will be too large!!!
     call matrix_pf_output(npfs-1,pf_mpi_comm,state_dim,cnt,pf%psi&
          &,pf%timestep,tsdata%is_analysis)
  end if
  
end subroutine output_from_pf

subroutine save_state(state,filename)
  use output_empire, only : unit_state,emp_e
  use sizes
  implicit none
  integer, parameter :: rk = kind(1.0d0)
  real(kind=rk), dimension(state_dim), intent(in) :: state
  character(256), intent(in) :: filename
  integer :: ios

  open(unit_state,file=trim(filename),iostat=ios,action='write',status='replace'&
       &,form='unformatted')
  if(ios .ne. 0)  then
     write(emp_e,*) 'PARTICLE FILTER DATA ERROR!!!!! Cannot open file '&
          &,filename
     write(emp_e,*) 'Very strange that I couldnt open it. Im going to stop&
          & now.'
     stop 'EMPIRE ERROR in SAVE_STATE'
  end if
  write(unit_state) state
  close(unit_state)
end subroutine save_state


subroutine get_state(state,filename)
  use output_empire, only : unit_state,emp_e
  use sizes
  implicit none
  integer, parameter :: rk = kind(1.0d0)
  real(kind=rk), dimension(state_dim), intent(out) :: state
  character(256), intent(in) :: filename
  integer :: ios

  open(unit_state,file=trim(filename),iostat=ios,action='read',status='old',form='un&
       &formatted')
  if(ios .ne. 0)  then
     write(emp_e,*) 'PARTICLE FILTER DATA ERROR!!!!! Cannot open file '&
          &,filename
     write(emp_e,*) 'Very strange that I couldnt open it. Im going to stop&
          & now.'
     stop
  end if
  read(unit_state) state
  close(unit_state)
end subroutine get_state