pf_control.f90

Go to the documentation of this file.

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! Time-stamp: <2016-11-23 11:47:49 pbrowne>
!!!
!!!    module to hold all the information to control the the main program
!!!    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
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

module pf_control
  implicit none
  type, public :: pf_control_type
     integer :: nens
     real(kind=kind(1.0D0)), allocatable, dimension(:) :: weight
     integer :: time_obs
     integer :: time_bwn_obs
     real(kind=kind(1.0D0)) :: nudgefac
     logical :: gen_data
     logical :: gen_q
     integer :: timestep=0
     real(kind=kind(1.0D0)), allocatable, dimension(:,:) :: psi
     real(kind=kind(1.0D0)), allocatable, dimension(:) :: mean
     real(kind=kind(1.0D0)) :: nfac
     real(kind=kind(1.0D0)) :: ufac
     real(kind=kind(1.0D0)) :: efac
     real(kind=kind(1.0D0)) :: keep
     real(kind=kind(1.0D0)) :: time
     real(kind=kind(1.0D0)) :: qscale
     real(kind=kind(1.0d0)) :: rho
     real(kind=kind(1.0d0)) :: len
     integer :: couple_root
     logical :: use_talagrand
     logical :: output_weights
     logical :: output_forecast
     logical :: use_mean
     logical :: use_variance
     logical :: use_traj
     logical :: use_spatial_rmse
     logical :: use_ens_rmse
     character(250) :: rmse_filename
     character(250) :: forecast_path
     character(50) :: relaxation_type
     real(kind=kind(1.0d0)) :: relaxation_freetime
     real(kind=kind(1.0d0)) :: power_law_p
     
     integer, dimension(:,:), allocatable :: talagrand
     integer :: count
     integer,allocatable, dimension(:) :: particles
     character(2) :: filter
     character(1) :: init
     
  end type pf_control_type
  type(pf_control_type), save :: pf

  
contains
  subroutine set_pf_controls
    use output_empire, only : emp_o
    
    write(emp_o,'(A)') 'Opening empire namelist file:'
    
    call parse_pf_parameters
    
    pf%efac = 0.001/pf%nens
    write(emp_o,'(A)') 'empire namelist successfully read to control pf code.'
    call flush(emp_o)
    
    
  end subroutine set_pf_controls



    subroutine parse_pf_parameters
      use var_data
      use output_empire, only : unit_nml,emp_e
      implicit none
      integer :: ios

      integer :: time_obs=-1
      integer :: time_bwn_obs=-1 
      real(kind=kind(1.0D0)) :: nudgefac=-1.0d0
      logical :: gen_data,gen_Q
      real(kind=kind(1.0D0)) :: nfac=-1.0d0                
      real(kind=kind(1.0D0)) :: ufac=-1.0d0
      real(kind=kind(1.0D0)) :: Qscale=-1.0d0
      real(kind=kind(1.0D0)) :: rho=0.0d0
      real(kind=kind(1.0d0)) :: len=-1.0d0
      real(kind=kind(1.0D0)) :: keep
      logical :: use_talagrand,use_mean,use_variance,use_traj&
           &,use_spatial_rmse,use_ens_rmse,output_weights,output_forecast
      character(2) :: filter='++'
      character(1) :: init='+'
      character(250) :: rmse_filename='rmse'
      character(250) :: forecast_path='forecast/'
      character(50) :: relaxation_type='zero_linear'
      real(kind=kind(1.0d0)) :: power_law_p=1.0d0
      real(kind=kind(1.0d0)) :: relaxation_freetime=0.6d0
      
      logical :: file_exists

      namelist/pf_params/time_obs,time_bwn_obs,&
      &nudgefac,& 
      &gen_data,gen_q,&
      &nfac,&
      &keep,&
      &ufac,&                
      &qscale,&
      &rho,&
      &len,&
      &use_talagrand,use_mean,use_variance,use_traj,use_spatial_rmse,use_ens_rmse,&
      &output_weights,&
      &output_forecast,&
      &filter,&
      &init,&
      &rmse_filename,&
      &forecast_path,&
      &relaxation_type,power_law_p,relaxation_freetime


      gen_data = .false.
      gen_q = .false.
      use_talagrand = .false.
      use_mean = .false.
      use_variance = .false.
      use_traj = .false.
      use_spatial_rmse = .false.
      use_ens_rmse = .false.
      output_weights=.false.
      output_forecast=.false.


      inquire(file='pf_parameters.dat',exist=file_exists)
      if(file_exists) then
         open(unit_nml,file='pf_parameters.dat',iostat=ios,action='read'&
              &,status='old')
         if(ios .ne. 0) then
            write(emp_e,*) 'Cannot open pf_parameters.dat'
            stop
         end if
      else
         inquire(file='empire.nml',exist=file_exists)
         if(file_exists) then
            open(unit_nml,file='empire.nml',iostat=ios,action='read'&
                 &,status='old')
            if(ios .ne. 0) stop 'Cannot open empire.nml'
         else
            write(emp_e,*) 'ERROR: cannot find pf_parameters.dat or empire.nml'
            stop '-1'
         end if
      end if
         
      read(unit_nml,nml=pf_params) 
      close(unit_nml)

      if(time_obs .gt. -1) then
         print*,'read time_obs = ',time_obs
         pf%time_obs = time_obs
      end if
      if(time_bwn_obs .gt. -1) then
         print*,'read time_bwn_obs = ',time_bwn_obs
         pf%time_bwn_obs = time_bwn_obs
      end if
      if(nudgefac .gt. -1.0d0) then
         print*,'read nudgefac = ',nudgefac
         pf%nudgefac = nudgefac
      end if

      if(nfac .gt. -1.0d0) then
         print*,'read nfac = ',nfac
         pf%nfac = nfac
      end if

      if(keep .gt. -1.0d0) then
         print*,'read keep = ',keep
         pf%keep = keep
      end if

      if(ufac .gt. -1.0d0) then
         print*,'read ufac = ',ufac
         pf%ufac = ufac
      end if

      !real(kind=kind(1.0D0)) :: ufac=-1.0d0
      if(qscale .gt. -1.0d0) then
         print*,'read Qscale = ',qscale
         pf%Qscale = qscale
      end if
      

      if(rho .gt. 0.0d0) then
         print*,'read rho = ',rho
         pf%rho = rho
      elseif(rho .lt. 0.0d0) then
         print*,'read rho = ',rho,' WARNING ABOUT THAT ONE! rho is nor&
              &mally positive'
         pf%rho = rho
      else
         pf%rho = rho
      end if


      if(len .ge. 0.0d0) then
         print*,'read len = ',len
         pf%len = len
       else
         pf%len = len
      end if

      if(relaxation_type .ne. 'zero_linear') then
         print*,'read relaxation_type = ',relaxation_type
         pf%relaxation_type = relaxation_type
      else
         pf%relaxation_type = relaxation_type
      end if

      if(relaxation_freetime .ne. 0.6d0) then
         print*,'read relaxation_freetime = ',relaxation_freetime
         pf%relaxation_freetime = relaxation_freetime
      else
         pf%relaxation_freetime = relaxation_freetime
      end if

      if(power_law_p .ne. 1.0d0) then
         print*,'read power_law_p = ',power_law_p
         pf%power_law_p = power_law_p
      else
         pf%power_law_p = power_law_p
      end if

      
      !logical ::
      !use_talagrand,use_mean,use_variance,use_traj,use_spatial_rmse
      
      if(use_spatial_rmse) then
         print*,'going to output spatial Root mean squared errors'
      end if

      if(use_ens_rmse) then
         print*,'going to output field of Root mean squared errors'
      end if
      
      if(use_mean) then
         print*,'going to output ensemble mean'
      end if

      if(use_variance) then
         print*,'going to output ensemble variance'
      end if
      
      if(output_weights) then
         print*,'going to output ensemble weights'
      end if

      if(output_forecast) then
         print*,'going to output ensemble forecasts'
      end if

      
      if(use_traj) then
         print*,'going to output trajectories'
      end if


      if(rmse_filename .ne. 'rmse') then
         print*,'read rmse_filename = ',rmse_filename
      end if
      pf%rmse_filename = rmse_filename

      if(forecast_path .ne. 'forecast/') then
         print*,'read forecast_path = ',forecast_path
      end if
      pf%forecast_path = forecast_path
      
      
      if(filter .ne. '++') then
         print*,'read filter = ',filter
         pf%filter = filter
      end if


      !ensure that if we are generating the data then SE is selected
      if(gen_data) then
         pf%filter = 'SE'         
      end if



      !let us verify pf%filter
      select case(pf%filter)
      case('EW')
         print*,'Running the equivalent weights particle filter'
      case('EZ')
         print*,'Running the Zhu equal weights particle filter'
      case('SE')
         print*,'Running a stochastic ensemble'
      case('DE')
         print*,'Running a deterministic ensemble'
      case('SI')
         print*,'Running the SIR particle filter'
      case('ET')
         print*,'filter read as ET. This is depreciated. Changing to L&
              &E'
         pf%filter = 'LE'
         print*,'Running the Local Ensemble Transform Kalman Filter'
         print*,'With random noise'
         print*,'For LETKF without random noise set filter="LD"'
      case('EA')
         print*,'Running the Ensemble Adjustment Kalman Filter'
         write(emp_e,*) 'Error: The EAKF is not implemented here yet'
         stop '-557'
      case('LE')
         print*,'Running the Local Ensemble Transform Kalman Filter'
         print*,'With random noise'
      case('LS')
         print*,'Running the Local Ensemble Transform Kalman Smoother'
         print*,'With random noise'
      case('LD')
         print*,'Running the Local Ensemble Transform Kalman Filter'
         print*,'With NO random noise'
      case('3D')
         print*,'Running a stochastic ensemble and 3DVar at observatio&
              &n times'
         call set_var_controls
      case default
         write(emp_e,*) 'Error: Incorrect filter type selected:', pf%filter
         write(emp_e,*) 'Please ensure that pf%filter in empire.nml is either:'
         write(emp_e,*) 'EW        the equivalent weights particle f&
              &ilter'
         write(emp_e,*) 'EZ        the Zhu equal weights particle f&
              &ilter'
         write(emp_e,*) 'SE        a stochastic ensemble'
         write(emp_e,*) 'DE        a deterministic ensemble'
         write(emp_e,*) 'SI        the SIR particle filter'
         write(emp_e,*) 'LE        the Local Ensemble Transform Kalman Filter'
         write(emp_e,*) '          with random noise'
         write(emp_e,*) 'LD        the Local Ensemble Transform Kalman Filter'
         write(emp_e,*) '          without random noise'
         write(emp_e,*) 'LS        the Local Ensemble Transform Kalman Smoother'
         write(emp_e,*) '          with random noise'
         write(emp_e,*) '3D        3DVar'
         stop
      end select

      
      if(init .ne. '+') then
         print*,'read init = ',init
         pf%init = init
      end if

      pf%gen_data = gen_data
      pf%gen_Q = gen_q
      pf%use_talagrand = use_talagrand
      pf%use_mean = use_mean
      pf%use_variance = use_variance
      pf%use_traj = use_traj
      pf%use_spatial_rmse = use_spatial_rmse
      pf%use_ens_rmse = use_ens_rmse
      pf%output_weights = output_weights
      pf%output_forecast = output_forecast
      
      !check for specific things we can't do if running the truth
      if(pf%gen_data) then
         pf%gen_Q = .false.
         pf%use_talagrand=.false.
         pf%use_spatial_rmse = .false.
         pf%use_ens_rmse = .false.
         pf%output_forecast = .false.
      end if

      
    end subroutine parse_pf_parameters


    
    subroutine deallocate_pf
      deallocate(pf%weight)
      deallocate(pf%psi)
      if(allocated(pf%talagrand)) deallocate(pf%talagrand)
      deallocate(pf%particles)
    end subroutine deallocate_pf


end module pf_control