v3fit_params.f

Go to the documentation of this file.

!*******************************************************************************
!
!  Note separating the Doxygen comment block here so detailed decription is
!  found in the Module not the file.
!
!*******************************************************************************
 
      MODULE v3fit_params
 
      USE stel_kinds, only: rprec
      USE data_parameters
      USE model
 
      IMPLICIT NONE
 
!*******************************************************************************
!  param module parameters
!*******************************************************************************
      INTEGER, PARAMETER :: param_range_no_type        = -1
      INTEGER, PARAMETER :: param_range_value_type     = 0
      INTEGER, PARAMETER :: param_range_infinity_type  = 1
      INTEGER, PARAMETER :: param_range_parameter_type = 2
 
      INTEGER, PARAMETER :: param_max_increment_steps = 5
      INTEGER, PARAMETER :: param_div_factor = 10.0
 
!*******************************************************************************
!  DERIVED-TYPE DECLARATIONS
!  1) parameter base class
!
!*******************************************************************************
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      TYPE param_recon_class
         REAL (rprec)                          :: vrnc = 0.0
 
!  Parameter ranges
         INTEGER, DIMENSION(2) :: range_type = param_range_no_type
         INTEGER, DIMENSION(2) :: range_id = data_no_id
         INTEGER, DIMENSION(2,data_max_indices) :: range_indices = 0
         REAL (rprec), DIMENSION(2) :: range_value = 0.0
 
         REAL (rprec)                          :: delta
 
         REAL (rprec), DIMENSION(:), POINTER   :: sem => null()
      END TYPE
 
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      TYPE param_locking_class
         INTEGER, DIMENSION(:), POINTER   :: ids => null()
         INTEGER, DIMENSION(:,:), POINTER :: indices => null()
      END TYPE
 
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      TYPE param_class
         INTEGER                              :: param_id = data_no_id
         INTEGER, DIMENSION(data_max_indices) :: indices = 0
 
         TYPE (param_recon_class), POINTER    :: recon => null()
         TYPE (param_locking_class), POINTER  :: locks => null()
 
         REAL (rprec)                         :: sigma = 0.0
         REAL (rprec), DIMENSION(:), POINTER  :: correlation => null()
      END TYPE
 
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      TYPE param_pointer
         TYPE (param_class), POINTER :: p => null()
      END TYPE
 
!*******************************************************************************
!  INTERFACE BLOCKS
!*******************************************************************************
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      INTERFACE param_construct
         MODULE PROCEDURE param_construct_basic,                               &
     &                    param_construct_recon,                               &
     &                    param_construct_locking
      END INTERFACE
 
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      INTERFACE param_write_step_data
         MODULE PROCEDURE param_write_step_data_1,                             &
     &                    param_write_step_data_2
      END INTERFACE
 
      CONTAINS
!*******************************************************************************
!  CONSTRUCTION SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      FUNCTION param_construct_basic(a_model, param_name, indices,             &
     &                               num_params)
 
      IMPLICIT NONE
 
!  Declare Arguments
      TYPE (param_class), POINTER       :: param_construct_basic
      TYPE (model_class), INTENT(in)    :: a_model
      CHARACTER (len=*), INTENT(in)     :: param_name
      INTEGER, DIMENSION(2), INTENT(in) :: indices
      INTEGER, INTENT(in)               :: num_params
 
!  local variables
      REAL (rprec)                      :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      ALLOCATE(param_construct_basic)
 
      param_construct_basic%param_id =                                         &
     &   model_get_param_id(a_model, trim(param_name))
 
      ALLOCATE(param_construct_basic%correlation(num_params))
      param_construct_basic%correlation = 0.0
 
      param_construct_basic%indices = indices
 
      CALL profiler_set_stop_time('param_construct_basic', start_time)
 
      END FUNCTION
 
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      FUNCTION param_construct_recon(a_model, param_name, indices,             &
     &                               vrnc, range_type, range_indices,          &
     &                               range_value, num_signals,                 &
     &                               num_params)
      USE v3_utilities
 
      IMPLICIT NONE
 
!  Declare Arguments
      TYPE (param_class), POINTER :: param_construct_recon
      TYPE (model_class), INTENT(in)              :: a_model
      CHARACTER (len=*), INTENT(in)               :: param_name
      INTEGER, DIMENSION(2), INTENT(in)           :: indices
      REAL (rprec), INTENT(in)                    :: vrnc
      CHARACTER (len=*), DIMENSION(2), INTENT(in) :: range_type
      INTEGER, DIMENSION(2,2), INTENT(in)         :: range_indices
      REAL (rprec), DIMENSION(2), INTENT(in)      :: range_value
      INTEGER, INTENT(in)                         :: num_signals
      INTEGER, INTENT(in)                         :: num_params
 
!  local variables
      REAL (rprec)                                :: value
      REAL (rprec)                                :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      param_construct_recon => param_construct(a_model, param_name,            &
     &                                         indices, num_params)
 
      IF (.not.model_is_recon_param(a_model,                                   &
     &            param_construct_recon%param_id)) THEN
         CALL err_fatal('param_construct: ' // trim(param_name) //             &
     &                  ' is an invalid reconstruction ' //                    &
     &                  'parameter')
      END IF
 
      ALLOCATE(param_construct_recon%recon)
      param_construct_recon%recon%vrnc = vrnc
 
      SELECT CASE (trim(range_type(1)))
 
         CASE ('value')
            param_construct_recon%recon%range_type(1) =                        &
     &         param_range_value_type
            param_construct_recon%recon%range_value(1) =                       &
     &         range_value(1)
        
         CASE ('infinity')
            param_construct_recon%recon%range_type(1) =                        &
     &         param_range_infinity_type
 
         CASE DEFAULT
            param_construct_recon%recon%range_type(1) =                        &
     &         param_range_parameter_type
            param_construct_recon%recon%range_id(1) =                          &
     &         model_get_param_id(a_model, trim(range_type(1)))
            param_construct_recon%recon%range_indices(1,:) =                   &
     &         range_indices(1,:)
 
      END SELECT
 
      SELECT CASE (trim(range_type(2)))
 
         CASE ('value')
            param_construct_recon%recon%range_type(2) =                        &
     &         param_range_value_type
            param_construct_recon%recon%range_value(2) =                       &
     &         range_value(2)
 
         CASE ('infinity')
            param_construct_recon%recon%range_type(2) =                        &
     &         param_range_infinity_type
 
         CASE DEFAULT
            param_construct_recon%recon%range_type(2) =                        &
     &         param_range_parameter_type
            param_construct_recon%recon%range_id(2) =                          &
     &         model_get_param_id(a_model, trim(range_type(2)))
            param_construct_recon%recon%range_indices(2,:) =                   &
     &         range_indices(2,:)
 
      END SELECT
 
      ALLOCATE(param_construct_recon%recon%sem(num_signals))
      param_construct_recon%recon%sem = 0.0
 
!  Check if the inital parameter value is in range.
      value = param_get_value(param_construct_recon, a_model)
      CALL assert(param_is_in_lower_range(param_construct_recon,               &
     &                                    a_model, value),                     &
     &            trim(param_get_name(param_construct_recon, a_model))         &
     &            // ' inital value is outside reconstruction parameter'       &
     &            //' lower bound.')
 
      CALL assert(param_is_in_upper_range(param_construct_recon,               &
     &                                    a_model, value),                     &
     &            trim(param_get_name(param_construct_recon, a_model))         &
     &            // ' inital value is outside reconstruction parameter'       &
     &            // ' upper bound.')
 
      CALL profiler_set_stop_time('param_construct_recon', start_time)
 
      END FUNCTION
 
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      FUNCTION param_construct_locking(a_model, param_name, indices,           &
     &                                 set, set_indices, set_coeff,            &
     &                                 eq_comm)
      USE v3_utilities
 
      IMPLICIT NONE
 
!  Declare Arguments
      TYPE (param_class), POINTER            :: param_construct_locking
      TYPE (model_class), INTENT(inout)      :: a_model
      CHARACTER (len=*), INTENT(in)          :: param_name
      INTEGER, DIMENSION(2), INTENT(in)      :: indices
      CHARACTER (len=*), DIMENSION(:), INTENT(in) :: set
      INTEGER, DIMENSION(:,:), INTENT(in)    :: set_indices
      REAL (rprec), DIMENSION(:), INTENT(in) :: set_coeff
      INTEGER, INTENT(in)                    :: eq_comm
 
!  local variables
      INTEGER                                :: i
      REAL (rprec)                           :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      param_construct_locking => param_construct(a_model, param_name,          &
     &                                           indices, SIZE(set))
 
      IF (.not.model_is_recon_param(a_model,                                   &
     &            param_construct_locking%param_id)) THEN
         CALL err_fatal('param_construct: ' // trim(param_name) //             &
     &                  ' is an invalid locking parameter')
      END IF
 
      ALLOCATE(param_construct_locking%locks)
 
      ALLOCATE(param_construct_locking%locks%ids(SIZE(set)))
      DO i = 1, SIZE(set)
         param_construct_locking%locks%ids(i) =                                &
     &      model_get_param_id(a_model, trim(set(i)))
      END DO
 
      ALLOCATE(param_construct_locking%locks%indices(SIZE(set),                &
     &                                               data_max_indices))
      param_construct_locking%locks%indices = set_indices
 
!  Use the correlation array as the array to hold the set coefficients. This
!  array was allocated by param_construct.
      param_construct_locking%correlation = set_coeff
 
!  Set the inital value of the locked parameter.
      CALL param_set_lock_value(param_construct_locking, a_model,              &
     &                          eq_comm)
 
      CALL profiler_set_stop_time('param_construct_locking', start_time)
 
      END FUNCTION
 
!*******************************************************************************
!  DESTRUCTION SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      SUBROUTINE param_destruct(this)
 
      IMPLICIT NONE
 
!  Declare Arguments
      TYPE (param_class), POINTER :: this
 
!  Start of executable code
      this%param_id = data_no_id
      this%indices = 0
 
      IF (ASSOCIATED(this%recon)) THEN
 
         IF (ASSOCIATED(this%recon%sem)) THEN
            DEALLOCATE(this%recon%sem)
            this%recon%sem => null()
         END IF
 
         DEALLOCATE(this%recon)
         this%recon => null()
      END IF
 
      IF (ASSOCIATED(this%locks)) THEN
 
         IF (ASSOCIATED(this%locks%ids)) THEN
            DEALLOCATE(this%locks%ids)
            this%locks%ids => null()
         END IF
 
         IF (ASSOCIATED(this%locks%indices)) THEN
            DEALLOCATE(this%locks%indices)
            this%locks%indices => null()
         END IF
 
         DEALLOCATE(this%locks)
         this%locks => null()
      END IF
 
      IF (ASSOCIATED(this%correlation)) THEN
         DEALLOCATE(this%correlation)
         this%correlation => null()
      END IF
 
      DEALLOCATE(this)
 
      END SUBROUTINE
 
!*******************************************************************************
!  SETTER SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      SUBROUTINE param_set_value(this, a_model, value, eq_comm,                &
     &                           is_central)
 
      IMPLICIT NONE
 
!  Declare Arguments
      TYPE (param_class), INTENT(in)    :: this
      TYPE (model_class), INTENT(inout) :: a_model
      REAL (rprec), INTENT(in)          :: value
      INTEGER, INTENT(in)               :: eq_comm
      LOGICAL, INTENT(in)               :: is_central
 
!  Local variables
      REAL (rprec)                      :: set_value
      REAL (rprec)                      :: off_set_value
      REAL (rprec)                      :: upper_value
      REAL (rprec)                      :: lower_value
      REAL (rprec)                      :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      set_value = value
 
      IF (is_central) THEN
         off_set_value = this%recon%vrnc/2.0
      ELSE
         off_set_value = 0.0
      END IF
 
!  Check if the value is within the lower range. If the value is outside,
!  truncate the value to the lower range.
      IF (this%recon%range_type(1) .ne. param_range_infinity_type) THEN
         lower_value = param_get_lower_range_value(this, a_model)
 
         IF (set_value - off_set_value .lt. lower_value) THEN
 
!  The set value must be at least the lower value plus the offset value. When
!  not using central differencing, the the offset is zero.
            set_value = lower_value + off_set_value
 
!  Need to check the upper range if using central differencing.
            IF (is_central .and. (this%recon%range_type(2) .ne.                &
     &                            param_range_infinity_type)) THEN
 
               upper_value = param_get_upper_range_value(this, a_model)
 
               IF (set_value + off_set_value .gt. upper_value) THEN
!  Both bounds are two narrow to fit the new value. Split the difference.
                  set_value = (upper_value + lower_value)/2.0
               END IF
            END IF
 
            CALL model_set_param(a_model, this%param_id,                       &
     &                           this%indices(1), this%indices(2),             &
     &                           set_value, eq_comm)
 
            CALL profiler_set_stop_time('profiler_get_start_time',             &
     &                                  start_time)
            RETURN
         END IF
      END IF
 
!  The value is in the lower range. Check upper range. If the value is outside,
!  truncate the value to the upper range.
      IF (this%recon%range_type(2) .ne. param_range_infinity_type) THEN
         upper_value = param_get_upper_range_value(this, a_model)
 
         IF (set_value + off_set_value .gt. upper_value) THEN
 
!  The set value must be no more than the upper value minus the offset value.
!  When not using central differencing, the the offset is zero. There is no need
!  to check the lower range because it has already been confirmed that the lower
!  range is negative infinity.
            set_value = upper_value - off_set_value
 
            CALL model_set_param(a_model, this%param_id,                       &
     &                           this%indices(1), this%indices(2),             &
     &                           set_value, eq_comm)
 
            CALL profiler_set_stop_time('profiler_get_start_time',             &
     &                                  start_time)
            RETURN
         END IF
      END IF
 
!  The value is in range. Set the value to this.
      CALL model_set_param(a_model, this%param_id,                             &
     &                     this%indices(1), this%indices(2),                   &
     &                     set_value, eq_comm)
 
      CALL profiler_set_stop_time('profiler_get_start_time', start_time)
 
      END SUBROUTINE
 
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      SUBROUTINE param_set_lock_value(this, a_model, eq_comm)
 
      IMPLICIT NONE
 
!  Declare Arguments
      TYPE (param_class), INTENT(in)    :: this
      TYPE (model_class), INTENT(inout) :: a_model
      INTEGER, INTENT(in)               :: eq_comm
 
!  Local variables
      INTEGER                           :: i
      REAL (rprec)                      :: temp
      REAL (rprec)                      :: inital_value
      REAL (rprec)                      :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      inital_value = model_get_param_value(a_model, this%param_id,             &
     &                                     this%indices(1),                    &
     &                                     this%indices(2))
 
      temp = 0.0
      DO i = 0, SIZE(this%locks%ids)
         temp = temp + this%correlation(i)                                     &
     &               * model_get_param_value(a_model,                          &
     &                                       this%locks%ids(i),                &
     &                                       this%locks%indices(i,1),          &
     &                                       this%locks%indices(i,2))
      END DO
 
!  Only update the lock value if that value changed.
      IF (inital_value .ne. temp) THEN
         CALL model_set_param(a_model, this%param_id,                          &
     &                        this%indices(1), this%indices(2),                &
     &                        temp, eq_comm)
      END iF
 
      CALL profiler_set_stop_time('param_set_lock_value', start_time)
 
      END SUBROUTINE
 
!*******************************************************************************
!  GETTER SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      FUNCTION param_get_value(this, a_model)
 
      IMPLICIT NONE
 
!  Declare Arguments
      REAL (rprec)                   :: param_get_value
      TYPE (param_class), INTENT(in) :: this
      TYPE (model_class), INTENT(in) :: a_model
 
!  local variables
      REAL (rprec)                   :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      param_get_value = model_get_param_value(a_model,                         &
     &                                        this%param_id,                   &
     &                                        this%indices(1),                 &
     &                                        this%indices(2))
 
      CALL profiler_set_stop_time('param_get_value', start_time)
 
      END FUNCTION
 
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      FUNCTION param_get_name(this, a_model)
 
      IMPLICIT NONE
 
!  Declare Arguments
      CHARACTER (len=data_name_length) :: param_get_name
      TYPE (param_class), INTENT(in)   :: this
      TYPE (model_class), INTENT(in)   :: a_model
 
!  local variables
      REAL (rprec)                     :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      param_get_name = model_get_param_name(a_model, this%param_id)
 
      CALL profiler_set_stop_time('param_get_name', start_time)
 
      END FUNCTION
 
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      FUNCTION param_get_lower_range_value(this, a_model)
 
      IMPLICIT NONE
 
!  Declare Arguments
      REAL (rprec)                   :: param_get_lower_range_value
      TYPE (param_class), INTENT(in) :: this
      TYPE (model_class), INTENT(in) :: a_model
 
!  local variables
      REAL (rprec)                   :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      SELECT CASE (this%recon%range_type(1))
 
         CASE DEFAULT
            param_get_lower_range_value = this%recon%range_value(1)
 
         CASE (param_range_parameter_type)
            param_get_lower_range_value =                                      &
     &         model_get_param_value(a_model, this%recon%range_id(1),          &
     &                               this%recon%range_indices(1,1),            &
     &                               this%recon%range_indices(1,2))
 
      END SELECT
 
      CALL profiler_set_stop_time('param_get_lower_range_value',               &
     &                            start_time)
 
      END FUNCTION
 
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      FUNCTION param_get_upper_range_value(this, a_model)
 
      IMPLICIT NONE
 
!  Declare Arguments
      REAL (rprec)                   :: param_get_upper_range_value
      TYPE (param_class), INTENT(in) :: this
      TYPE (model_class), INTENT(in) :: a_model
 
!  local variables
      REAL (rprec)                   :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      SELECT CASE (this%recon%range_type(2))
 
         CASE DEFAULT
            param_get_upper_range_value = this%recon%range_value(2)
 
         CASE (param_range_parameter_type)
            param_get_upper_range_value =                                      &
     &         model_get_param_value(a_model, this%recon%range_id(2),          &
     &                               this%recon%range_indices(2,1),            &
     &                               this%recon%range_indices(2,2))
 
      END SELECT
 
      CALL profiler_set_stop_time('param_get_upper_range_value',               &
     &                            start_time)
 
      END FUNCTION
 
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      FUNCTION param_get_lower_range_type(this, a_model)
 
      IMPLICIT NONE
 
!  Declare Arguments
      CHARACTER (len=data_name_length) :: param_get_lower_range_type
      TYPE (param_class), INTENT(in)   :: this
      TYPE (model_class), INTENT(in)   :: a_model
 
!  local variables
      REAL (rprec)                     :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      SELECT CASE (this%recon%range_type(1))
 
         CASE DEFAULT
            param_get_lower_range_type = 'No Type'
 
         CASE (param_range_value_type)
            param_get_lower_range_type = 'value'
 
         CASE (param_range_infinity_type)
            param_get_lower_range_type = 'infinity'
 
         CASE (param_range_parameter_type)
            param_get_lower_range_type =                                       &
     &         model_get_param_name(a_model, this%recon%range_id(1))
 
      END SELECT
 
      CALL profiler_set_stop_time('param_get_lower_range_type',                &
     &                            start_time)
 
      END FUNCTION
 
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      FUNCTION param_get_upper_range_type(this, a_model)
 
      IMPLICIT NONE
 
!  Declare Arguments
      CHARACTER (len=data_name_length) :: param_get_upper_range_type
      TYPE (param_class), INTENT(in)   :: this
      TYPE (model_class), INTENT(in)   :: a_model
 
!  local variables
      REAL (rprec)                     :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      SELECT CASE (this%recon%range_type(2))
 
         CASE DEFAULT
            param_get_upper_range_type = 'No Type'
 
         CASE (param_range_value_type)
            param_get_upper_range_type = 'value'
 
         CASE (param_range_infinity_type)
            param_get_upper_range_type = 'infinity'
 
         CASE (param_range_parameter_type)
            param_get_upper_range_type =                                       &
     &         model_get_param_name(a_model, this%recon%range_id(2))
 
      END SELECT
 
      CALL profiler_set_stop_time('param_get_upper_range_type',                &
     &                            start_time)
 
      END FUNCTION
 
!*******************************************************************************
!  QUERY SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      FUNCTION param_is_in_lower_range(this, a_model, value)
 
      IMPLICIT NONE
 
!  Declare Arguments
      LOGICAL                        :: param_is_in_lower_range
      TYPE (param_class), INTENT(in) :: this
      TYPE (model_class), INTENT(in) :: a_model
      REAL (rprec), INTENT(in)       :: value
 
!  local variables
      REAL (rprec)                   :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      param_is_in_lower_range = .true.
      IF (this%recon%range_type(1) .ne. param_range_infinity_type) THEN
         param_is_in_lower_range =                                             &
     &      value .ge. param_get_lower_range_value(this, a_model)
      END IF
 
      CALL profiler_set_stop_time('param_is_in_lower_range', start_time)
 
      END FUNCTION
 
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      FUNCTION param_is_in_upper_range(this, a_model, value)
 
      IMPLICIT NONE
 
!  Declare Arguments
      LOGICAL                        :: param_is_in_upper_range
      TYPE (param_class), INTENT(in) :: this
      TYPE (model_class), INTENT(in) :: a_model
      REAL (rprec), INTENT(in)       :: value
 
!  local variables
      REAL (rprec)                   :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      param_is_in_upper_range = .true.
      IF (this%recon%range_type(2) .ne. param_range_infinity_type) THEN
         param_is_in_upper_range =                                             &
     &      value .le. param_get_upper_range_value(this, a_model)
      END IF
 
      CALL profiler_set_stop_time('param_is_in_upper_range', start_time)
 
      END FUNCTION
 
!*******************************************************************************
!  UTILITY SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      SUBROUTINE param_increment(this, a_model, eq_comm, is_central)
      USE v3_utilities, only: err_fatal
 
      IMPLICIT NONE
 
!  Declare Arguments
      TYPE (param_class), INTENT(inout) :: this
      TYPE (model_class), INTENT(inout) :: a_model
      INTEGER, INTENT(in)               :: eq_comm
      LOGICAL, INTENT(in)               :: is_central
 
!  local variables
      REAL (rprec)                      :: new_value
      REAL (rprec)                      :: value
      REAL (rprec)                      :: step
      INTEGER                           :: itry
      REAL (rprec)                      :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      value = param_get_value(this, a_model)
 
!  Take the absolute value here to control the direction of the step.
      IF (is_central) THEN
         step = abs(this%recon%vrnc)/2.0
      ELSE
         step = abs(this%recon%vrnc)
      END IF
 
      DO itry = 1, param_max_increment_steps
 
!  Test to make sure there is a step
         IF (step .eq. 0.0) THEN
            CALL err_fatal('param_increment: tried to take a ' //              &
     &                     'step size of 0')
         END IF
 
         IF (is_central) THEN
!  Test both upper and lower range.
            IF (param_is_in_upper_range(this, a_model,                         &
     &                                  value + step) .and.                    &
     &          param_is_in_lower_range(this, a_model,                         &
     &                                  value - step)) THEN
               new_value = value + step
 
!  Store the actual change in value.
               this%recon%delta = step
               EXIT
            END IF
 
         ELSE
 
!  Try a step in the positive direction. If the parameter is still in range set
!  otherwise.
            new_value = value + step
            IF (param_is_in_upper_range(this, a_model, new_value)) THEN
!  Store the actual change in value.
               this%recon%delta = step
               EXIT
            END IF
 
!  Try a step in the negative direction.
            new_value = value - step
            IF (param_is_in_lower_range(this, a_model, new_value)) THEN
!  Store the actual change in value.
               this%recon%delta = -step
               EXIT
            END IF
         END IF
 
!  Decrease the step and try again.
         step = step/param_div_factor
      END DO
 
      CALL model_set_param(a_model, this%param_id, this%indices(1),            &
     &                     this%indices(2), new_value, eq_comm)
 
      IF (itry .gt. param_max_increment_steps) THEN
         CALL err_fatal('param_increment: failed to change ' //                &
     &                  'reconstruction. Try decreaing rp_vrnc ' //            &
     &                  'or expanding the range.')
      END IF
 
      CALL profiler_set_stop_time('param_increment', start_time)
 
      END SUBROUTINE
 
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      SUBROUTINE param_decrement(this, a_model, eq_comm)
 
      IMPLICIT NONE
 
!  Declare Arguments
      TYPE (param_class), INTENT(inout) :: this
      TYPE (model_class), INTENT(inout) :: a_model
      INTEGER, INTENT(in)               :: eq_comm
 
!  local variables
      REAL (rprec)                      :: value
      REAL (rprec)                      :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      value = param_get_value(this, a_model)
      value = value - this%recon%delta
 
      CALL model_set_param(a_model, this%param_id, this%indices(1),            &
     &                     this%indices(2), value, eq_comm)
 
      CALL profiler_set_stop_time('param_decrement', start_time)
 
      END SUBROUTINE
 
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      SUBROUTINE param_write(this, iou, index, a_model)
 
      IMPLICIT NONE
 
!  Declare Arguments
      TYPE (param_class), INTENT(in) :: this
      INTEGER, INTENT(in)            :: iou
      INTEGER, INTENT(in)            :: index
      TYPE (model_class), INTENT(in) :: a_model
 
!  local variables
      REAL (rprec)                   :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      WRITE (iou,1000) index,                                                  &
     &                 param_get_name(this, a_model),                          &
     &                 this%indices(1), this%indices(2),                       &
     &                 param_get_value(this, a_model),                         &
     &                 this%sigma, this%recon%vrnc,                            &
     &                 param_get_upper_range_type(this, a_model),              &
     &                 param_get_lower_range_type(this, a_model),              &
     &                 param_get_upper_range_value(this, a_model),             &
     &                 param_get_lower_range_value(this, a_model),             &
     &                 this%recon%range_indices(2,:),                          &
     &                 this%recon%range_indices(1,:)
1000  FORMAT(i3,2x,a18,2(1x,i4),3(2x,es12.5),2(2x,a13),2(2x,es12.5),           &
     &       2(2x,i8,1x,i8))
 
      CALL profiler_set_stop_time('param_write', start_time)
 
      END SUBROUTINE
 
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      SUBROUTINE param_write_short(this, iou, index, a_model)
 
      IMPLICIT NONE
 
!  Declare Arguments
      TYPE (param_class), INTENT(in) :: this
      INTEGER, INTENT(in)            :: iou
      INTEGER, INTENT(in)            :: index
      TYPE (model_class), INTENT(in) :: a_model
 
!  local variables
      REAL (rprec)                   :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      WRITE (iou,1000) index,                                                  &
     &                 param_get_name(this, a_model),                          &
     &                 this%indices(1), this%indices(2),                       &
     &                 param_get_value(this, a_model), this%sigma
1000  FORMAT(i3,2x,a18,2(1x,i4),2(2x,es12.5))
 
      CALL profiler_set_stop_time('param_write_short', start_time)
 
      END SUBROUTINE
 
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      SUBROUTINE param_write_header(iou)
 
      IMPLICIT NONE
 
!  Declare Arguments
      INTEGER, INTENT(in) :: iou
 
!  local variables
      REAL (rprec)        :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      WRITE (iou,*)
      WRITE (iou,*) ' *** Reconstruction parameters'
      WRITE (iou,1000)
 
      CALL profiler_set_stop_time('param_write_header', start_time)
 
1000  FORMAT ('irp',2x,'p_type',13x,'inx1',1x,'inx2',2x,'value',9x,            &
     &        'sigma',9x,'vrnc',10x,'range_type_h',3x,'range_type_l',3x,       &
     &        'r_value_h',5x,'r_value_l',5x,                                   &
     &        'r_inx1_h',1x,'r_inx2_h',2x,'r_inx1_l',1x,'r_inx2_l')
 
      END SUBROUTINE
 
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      SUBROUTINE param_write_header_short(iou)
 
      IMPLICIT NONE
 
!  Declare Arguments
      INTEGER, INTENT(in) :: iou
 
!  local variables
      REAL (rprec)        :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      WRITE (iou,*)
      WRITE (iou,*) ' *** Derived parameters'
      WRITE (iou,1000)
 
      CALL profiler_set_stop_time('param_write_header_short',                  &
     &                            start_time)
 
1000  FORMAT ('irp',2x,'p_type',13x,'inx1',1x,'inx2',2x,'value',9x,            &
     &        'sigma')
 
      END SUBROUTINE
 
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      SUBROUTINE param_write_correlation(this, iou, a_model)
 
      IMPLICIT NONE
 
!  Declare Arguments
      TYPE (param_class), INTENT(in) :: this
      INTEGER, INTENT(in)            :: iou
      TYPE (model_class), INTENT(in) :: a_model
 
!  local variables
      CHARACTER (len=20)             :: row_format
      REAL (rprec)                   :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      WRITE (row_format,1000) SIZE(this%correlation)
      WRITE (iou,row_format) param_get_name(this, a_model),                    &
     &                       this%correlation
 
      CALL profiler_set_stop_time('param_write_correlation', start_time)
 
1000  FORMAT ('(a12,',i3,'(2x,es12.5))')
 
      END SUBROUTINE
 
!*******************************************************************************
!  NETCDF SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      SUBROUTINE param_write_step_data_1(this, a_model, result_ncid,           &
     &                                   current_step, index,                  &
     &                                   param_value_id, param_sigma_id,       &
     &                                   param_corr_id, param_sem_id)
      USE ezcdf
 
      IMPLICIT NONE
 
!  Declare Arguments
      TYPE (param_class), INTENT(in) :: this
      TYPE (model_class), INTENT(in) :: a_model
      INTEGER, INTENT(in)            :: result_ncid
      INTEGER, INTENT(in)            :: current_step
      INTEGER, INTENT(in)            :: index
      INTEGER, INTENT(in)            :: param_value_id
      INTEGER, INTENT(in)            :: param_sigma_id
      INTEGER, INTENT(in)            :: param_corr_id
      INTEGER, INTENT(in)            :: param_sem_id
 
!  local variables
      INTEGER                        :: status
      REAL (rprec)                   :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      CALL param_write_step_data(this, a_model, result_ncid,                   &
     &                           current_step, index,                          &
     &                           param_value_id, param_sigma_id,               &
     &                           param_corr_id)
 
      status = nf_put_vara_double(result_ncid, param_sem_id,                   &
     &                            (/ 1, index, current_step /),                &
     &                            (/ SIZE(this%recon%sem), 1, 1 /),            &
     &                            this%recon%sem)
 
      CALL profiler_set_stop_time('param_write_step_data_1', start_time)
 
      END SUBROUTINE
 
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      SUBROUTINE param_write_step_data_2(this, a_model, result_ncid,           &
     &                                   current_step, index,                  &
     &                                   param_value_id, param_sigma_id,       &
     &                                   param_corr_id)
      USE ezcdf
 
      IMPLICIT NONE
 
!  Declare Arguments
      TYPE (param_class), INTENT(in) :: this
      TYPE (model_class), INTENT(in) :: a_model
      INTEGER, INTENT(in)            :: result_ncid
      INTEGER, INTENT(in)            :: current_step
      INTEGER, INTENT(in)            :: index
      INTEGER, INTENT(in)            :: param_value_id
      INTEGER, INTENT(in)            :: param_sigma_id
      INTEGER, INTENT(in)            :: param_corr_id
 
!  local variables
      INTEGER                        :: status
      REAL (rprec)                   :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      status = nf_put_var1_double(result_ncid, param_value_id,                 &
     &                            (/ index, current_step /),                   &
     &                            param_get_value(this, a_model))
      CALL assert_eq(status, nf_noerr, nf_strerror(status))
 
      status = nf_put_var1_double(result_ncid, param_sigma_id,                 &
     &                            (/ index, current_step /),                   &
     &                            this%sigma)
      CALL assert_eq(status, nf_noerr, nf_strerror(status))
 
      status = nf_put_vara_double(result_ncid, param_corr_id,                  &
     &                            (/ 1, index, current_step /),                &
     &                            (/ SIZE(this%correlation), 1, 1 /),          &
     &                            this%correlation)
      CALL assert_eq(status, nf_noerr, nf_strerror(status))
 
      CALL profiler_set_stop_time('param_write_step_data_2', start_time)
 
      END SUBROUTINE
 
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      SUBROUTINE param_restart(this, a_model, result_ncid, current_step,       &
     &                         index, param_value_id, param_sigma_id,          &
     &                         param_corr_id, eq_comm, is_central)
      USE ezcdf
 
      IMPLICIT NONE
 
!  Declare Arguments
      TYPE (param_class), INTENT(inout) :: this
      TYPE (model_class), INTENT(inout) :: a_model
      INTEGER, INTENT(in)               :: result_ncid
      INTEGER, INTENT(in)               :: current_step
      INTEGER, INTENT(in)               :: index
      INTEGER, INTENT(in)               :: param_value_id
      INTEGER, INTENT(in)               :: param_sigma_id
      INTEGER, INTENT(in)               :: param_corr_id
      INTEGER, INTENT(in)               :: eq_comm
      LOGICAL, INTENT(in)               :: is_central
 
!  local variables
      INTEGER                           :: status
      REAL (rprec)                      :: temp_value
      REAL (rprec)                      :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      status = nf_get_var1_double(result_ncid, param_value_id,                 &
     &                            (/ index, current_step /), temp_value)
      CALL assert_eq(status, nf_noerr, nf_strerror(status))
      CALL param_set_value(this, a_model, temp_value, eq_comm,                 &
     &                     is_central)
 
      status = nf_put_var1_double(result_ncid, param_sigma_id,                 &
     &                            (/ index, current_step /),                   &
     &                            this%sigma)
      CALL assert_eq(status, nf_noerr, nf_strerror(status))
 
      status = nf_put_vara_double(result_ncid, param_corr_id,                  &
     &                            (/ 1, index, current_step /),                &
     &                            (/ SIZE(this%correlation), 1, 1 /),          &
     &                            this%correlation)
      CALL assert_eq(status, nf_noerr, nf_strerror(status))
 
      CALL profiler_set_stop_time('param_restart', start_time)
 
      END SUBROUTINE
 
!*******************************************************************************
!  MPI SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      SUBROUTINE param_sync_value(this, a_model, recon_comm, eq_comm,          &
     &                            is_central)
 
      IMPLICIT NONE
 
!  Declare Arguments
      TYPE (param_class), INTENT(inout) :: this
      TYPE (model_class), INTENT(inout) :: a_model
      INTEGER, INTENT(in)               :: recon_comm
      INTEGER, INTENT(in)               :: eq_comm
      LOGICAL, INTENT(in)               :: is_central
 
#if defined(MPI_OPT)
!  local variables
      INTEGER                           :: error
      REAL (rprec)                      :: value
      INTEGER                           :: mpi_rank
      REAL (rprec)                      :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      CALL mpi_comm_rank(recon_comm, mpi_rank, error)
 
      IF (mpi_rank .eq. 0) THEN
         value = param_get_value(this, a_model)
      END IF
 
      CALL mpi_bcast(value, 1, mpi_real8, 0, recon_comm, error)
 
      IF (mpi_rank .gt. 0) THEN
         CALL param_set_value(this, a_model, value, eq_comm, is_central)
      END IF
 
      CALL profiler_set_stop_time('param_sync_value', start_time)
 
#endif
 
      END SUBROUTINE
 
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      SUBROUTINE param_send_delta(this, index, recon_comm)
 
      IMPLICIT NONE
 
!  Declare Arguments
      TYPE (param_class), INTENT(inout) :: this
      INTEGER, INTENT(in)               :: index
      INTEGER, INTENT(in)               :: recon_comm
 
#if defined(MPI_OPT)
!  local variables
      INTEGER                           :: error
      INTEGER                           :: mpi_size
      INTEGER                           :: mpi_request
      REAL (rprec)                      :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      CALL mpi_comm_size(recon_comm, mpi_size, error)
 
!  Must use a non blocking send here because we are sending messages from the
!  main process to itself or need to send multiple sends before the first
!  recieve. Other wise it could block until a recieve request is called.
      CALL mpi_isend(this%recon%delta, 1, mpi_real8, 0, index,                 &
     &               recon_comm, mpi_request, error)
 
      CALL profiler_set_stop_time('param_send_delta', start_time)
#endif
 
      END SUBROUTINE
 
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      SUBROUTINE param_recv_delta(this, index, recon_comm)
 
      IMPLICIT NONE
 
!  Declare Arguments
      TYPE (param_class), INTENT(inout) :: this
      INTEGER, INTENT(in)               :: index
      INTEGER, INTENT(in)               :: recon_comm
 
#if defined(MPI_OPT)
!  local variables
      INTEGER                           :: error
      REAL (rprec)                      :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      CALL mpi_recv(this%recon%delta, 1, mpi_real8, mpi_any_source,            &
     &              index, recon_comm, mpi_status_ignore, error)
 
      CALL profiler_set_stop_time('param_recv_delta', start_time)
 
#endif
 
      END SUBROUTINE
 
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      SUBROUTINE param_sync_delta(this, recon_comm)
 
      IMPLICIT NONE
 
!  Declare Arguments
      TYPE (param_class), INTENT(inout) :: this
      INTEGER, INTENT(in)               :: recon_comm
 
#if defined(MPI_OPT)
!  local variables
      INTEGER                           :: error
      REAL (rprec)                      :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      CALL mpi_bcast(this%recon%delta, 1, mpi_real8, 0, recon_comm,            &
     &               error)
 
      CALL profiler_set_stop_time('param_sync_delta', start_time)
 
#endif
 
      END SUBROUTINE
 
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      SUBROUTINE param_sync_child(this, a_model, index, recon_comm,            &
     &                            eq_comm, is_central)
 
      IMPLICIT NONE
 
!  Declare Arguments
      TYPE (param_class), INTENT(inout) :: this
      TYPE (model_class), INTENT(inout) :: a_model
      INTEGER, INTENT(in)               :: index
      INTEGER, INTENT(in)               :: recon_comm
      INTEGER, INTENT(in)               :: eq_comm
      LOGICAL, INTENT(in)               :: is_central
 
#if defined(MPI_OPT)
!  local variables
      INTEGER                           :: error
      REAL (rprec)                      :: value
      INTEGER                           :: mpi_rank
      REAL (rprec)                      :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      CALL mpi_comm_rank(recon_comm, mpi_rank, error)
 
      IF (mpi_rank .eq. index) THEN
         value = param_get_value(this, a_model)
         CALL mpi_ssend(value, 1, mpi_real8, 0, mpi_rank, recon_comm,          &
     &                  error)
      ELSE IF (mpi_rank .eq. 0) THEN
         CALL mpi_recv(value, 1, mpi_real8, index, index, recon_comm,          &
     &                 mpi_status_ignore, error)
         CALL param_set_value(this, a_model, value, eq_comm, is_central)
      END IF
 
      CALL profiler_set_stop_time('param_sync_child', start_time)
 
#endif
 
      END SUBROUTINE
 
      END MODULE