mod_da_param_estimation.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scalelib.h"
module mod_da_param_estimation
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prof
 
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: da_param_estimation_setup
  public :: da_param_estimation_finalize
  public :: da_param_estimation_update
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  private :: estimation_getvar
  private :: estimation_putvar
  private :: estimation_log
  private :: phys2func
  private :: func2phys
 
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  integer, parameter :: MAXNUM_PARAM_ESTIMATION = 100
 
  integer :: NUM_PARAM_ESTIMATION
  logical :: DO_PARAM_ESTIMATION
 
  character(len=H_SHORT) :: ESTIMATION_TARGET(MAXNUM_PARAM_ESTIMATION)
 
  logical :: PEST_TRANSLATION(MAXNUM_PARAM_ESTIMATION)
 
  real(RP) :: PEST_ULIMIT(MAXNUM_PARAM_ESTIMATION)
  real(RP) :: PEST_LLIMIT(MAXNUM_PARAM_ESTIMATION)
 
  real(RP), allocatable :: param_array(:,:)
 
  integer :: datatype
 
  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine da_param_estimation_setup
    use mpi
    use scale_prc, only: &
       prc_abort
    use scale_comm_ensemble, only: &
       ensemble_nprocs => comm_ensemble_nprocs, &
       ensemble_myrank => comm_ensemble_myrank
    implicit none
 
    namelist / param_da_param_estimation / &
       estimation_target, &
       pest_translation,  &
       pest_ulimit,       &
       pest_llimit
 
    integer :: v
    integer :: ierr
    !---------------------------------------------------------------------------
 
    num_param_estimation = 0
    do_param_estimation = .false.
 
    estimation_target(:) = ''
 
    pest_translation(:) = .false.
 
    pest_ulimit(:) = huge(0.0_rp)
    pest_llimit(:) = -huge(0.0_rp)
 
    !--- gather parameters from ensemble communication
    if     ( rp == sp ) then
       datatype = mpi_real
    else if( rp == dp ) then
       datatype = mpi_double_precision
    else
       log_error("DA_param_estimation_setup",*) 'The precision has not been implemented yet:', rp
       call prc_abort
    endif
 
    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_da_param_estimation,iostat=ierr)
    if( ierr < 0 ) then !--- missing
       log_info("DA_param_estimation_setup",*) 'Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       log_error("DA_param_estimation_setup",*) 'Not appropriate names in namelist PARAM_DA_PARAM_ESTIMATION. Check!'
       call prc_abort
    endif
    log_nml(param_da_param_estimation)
 
    do v = 1, maxnum_param_estimation
      if( trim(estimation_target(v)) /= '' ) then
        num_param_estimation = num_param_estimation + 1
      end if
    end do
 
    if( num_param_estimation > maxnum_param_estimation ) then
       log_error("DA_param_estimation_setup",*) 'NUM_PARAM_ESTIMATION has exceeded the limit! ( MAX = ', maxnum_param_estimation, ')'
       call prc_abort
    end if
 
    if( num_param_estimation > 0 ) then
      log_info("DA_param_estimation_setup",*) 'Set up Parameter Estimation module. The number of estimation target:', num_param_estimation
      do_param_estimation = .true.
 
      allocate( param_array( ensemble_nprocs, num_param_estimation ) )
 
      call estimation_getvar
      call estimation_log
 
    end if
 
    return
  end subroutine da_param_estimation_setup
 
  !-----------------------------------------------------------------------------
  subroutine da_param_estimation_finalize
    implicit none
    !---------------------------------------------------------------------------
 
    if( allocated( param_array ) ) deallocate( param_array )
 
    return
  end subroutine da_param_estimation_finalize
 
  !-----------------------------------------------------------------------------
  subroutine da_param_estimation_update
    use scale_comm_ensemble, only: &
      ensemble_nprocs => comm_ensemble_nprocs
    use scale_letkf, only: &
      letkf_param_estimation_system
    implicit none
 
    integer :: n, v
    !---------------------------------------------------------------------------
 
    if( .NOT. do_param_estimation ) return
 
    !--- get parameters for estimation
    call estimation_getvar
 
    do v = 1, num_param_estimation
    do n = 1, ensemble_nprocs
       if( pest_translation(v) ) then ! encoding function by Kotsuki (2015)
          param_array(n,v) = phys2func( v, param_array(n,v) )
       end if
    end do
    end do
 
    !--- execute parameter estimation
    call letkf_param_estimation_system( num_param_estimation, &
                                        param_array(:,:)      )
 
    do v = 1, num_param_estimation
    do n = 1, ensemble_nprocs
       if( pest_translation(v) ) then ! decoding function by Kotsuki (2015)
          param_array(n,v) = func2phys( v, param_array(n,v) )
       end if
    end do
    end do
 
    !--- overwrite parameters
    call estimation_putvar
    call estimation_log
 
    return
  end subroutine da_param_estimation_update
 
  subroutine estimation_getvar
    use mpi
    use scale_prc, only: &
      prc_abort
    use scale_const, only: &
      undef => const_undef
    use scale_comm_ensemble, only: &
      ensemble_world  => comm_ensemble_world,  &
      ensemble_nprocs => comm_ensemble_nprocs, &
      ensemble_myrank => comm_ensemble_myrank
    use scale_atmos_phy_mp_tomita08, only: &
      atmos_phy_mp_tomita08_getvar
    use scale_atmos_phy_cp_kf, only: &
      atmos_phy_cp_kf_getvar
    implicit none
 
    real(rp), allocatable :: send(:)
    real(rp), allocatable :: recv(:)
 
    integer :: n, v
    integer :: ierr
    !---------------------------------------------------------------------------
 
    do v = 1, num_param_estimation
      select case( trim(estimation_target(v)) )
      case('MP_tomita08_drag_g')
         call atmos_phy_mp_tomita08_getvar( out_drag_g=param_array(1+ensemble_myrank,v) )
      case('MP_tomita08_re_qc')
         call atmos_phy_mp_tomita08_getvar( out_re_qc=param_array(1+ensemble_myrank,v) )
      case('MP_tomita08_re_qi')
         call atmos_phy_mp_tomita08_getvar( out_re_qi=param_array(1+ensemble_myrank,v) )
      case('MP_tomita08_Cr')
         call atmos_phy_mp_tomita08_getvar( out_cr=param_array(1+ensemble_myrank,v) )
      case('MP_tomita08_Cs')
         call atmos_phy_mp_tomita08_getvar( out_cs=param_array(1+ensemble_myrank,v) )
      case('CP_kf_delcape')
         call atmos_phy_cp_kf_getvar( out_delcape=param_array(1+ensemble_myrank,v) )
      case('CP_kf_deeplifetime')
         call atmos_phy_cp_kf_getvar( out_deeplifetime=param_array(1+ensemble_myrank,v) )
      case('CP_kf_shallowlifetime')
         call atmos_phy_cp_kf_getvar( out_shallowlifetime=param_array(1+ensemble_myrank,v) )
      case default
         log_error("DA_param_estimation_update",*) 'Invalid ESTIMATION_TARGET:', trim(estimation_target(v))
         call prc_abort
      end select
    end do
 
    allocate( send(num_param_estimation*ensemble_nprocs) )
    allocate( recv(num_param_estimation*ensemble_nprocs) )
 
    send(:) = undef
    recv(:) = undef
 
    do v = 1, num_param_estimation
       send( v + ensemble_myrank*num_param_estimation ) = param_array(1+ensemble_myrank,v)
    end do
 
    n = ensemble_myrank * num_param_estimation + 1
 
    call mpi_allgather( send(n),              &
                        num_param_estimation, &
                        datatype,             &
                        recv(1),              &
                        num_param_estimation, &
                        datatype,             &
                        ensemble_world,       &
                        ierr                  )
 
    do v = 1, num_param_estimation
    do n = 1, ensemble_nprocs
       param_array(n,v) = recv( v + (n-1)*num_param_estimation )
    end do
    end do
 
    deallocate( send )
    deallocate( recv )
 
    return
  end subroutine estimation_getvar
 
  subroutine estimation_putvar
    use scale_prc, only: &
      prc_abort
    use scale_comm_ensemble, only: &
      ensemble_myrank => comm_ensemble_myrank
    use scale_atmos_phy_mp_tomita08, only: &
      atmos_phy_mp_tomita08_putvar
    use scale_atmos_phy_cp_kf, only: &
      atmos_phy_cp_kf_putvar
    implicit none
 
    integer :: v
    !---------------------------------------------------------------------------
 
    do v = 1, num_param_estimation
      select case( trim(estimation_target(v)) )
      case('MP_tomita08_drag_g')
         call atmos_phy_mp_tomita08_putvar( in_drag_g=param_array(1+ensemble_myrank,v) )
      case('MP_tomita08_re_qc')
         call atmos_phy_mp_tomita08_putvar( in_re_qc=param_array(1+ensemble_myrank,v) )
      case('MP_tomita08_re_qi')
         call atmos_phy_mp_tomita08_putvar( in_re_qi=param_array(1+ensemble_myrank,v) )
      case('MP_tomita08_Cr')
         call atmos_phy_mp_tomita08_putvar( in_cr=param_array(1+ensemble_myrank,v) )
      case('MP_tomita08_Cs')
         call atmos_phy_mp_tomita08_putvar( in_cs=param_array(1+ensemble_myrank,v) )
      case('CP_kf_delcape')
         call atmos_phy_cp_kf_putvar( in_delcape=param_array(1+ensemble_myrank,v) )
      case('CP_kf_deeplifetime')
         call atmos_phy_cp_kf_putvar( in_deeplifetime=param_array(1+ensemble_myrank,v) )
      case('CP_kf_shallowlifetime')
         call atmos_phy_cp_kf_putvar( in_shallowlifetime=param_array(1+ensemble_myrank,v) )
      case default
         log_error("DA_param_estimation_update",*) 'Invalid ESTIMATION_TARGET:', trim(estimation_target(v))
         call prc_abort
      end select
    end do
 
    return
  end subroutine estimation_putvar
 
  subroutine estimation_log
    use scale_const, only: &
      undef => const_undef
    use scale_comm_ensemble, only: &
      ensemble_nprocs => comm_ensemble_nprocs, &
      ensemble_myrank => comm_ensemble_myrank
    use scale_time, only: &
      time_nowdate
    use scale_statistics, only: &
      average => statistics_average, &
      stddev  => statistics_stddev
    implicit none
 
    character(len=20) :: nowdate
 
    integer :: i, v
    !---------------------------------------------------------------------------
 
    do v = 1, num_param_estimation
      nowdate = '****/**/**-**:**:** '
      write( nowdate( 1: 4), '(i4)' ) time_nowdate(1)
      write( nowdate( 6: 7), '(i2)' ) time_nowdate(2)
      write( nowdate( 9:10), '(i2)' ) time_nowdate(3)
      write( nowdate(12:13), '(i2)' ) time_nowdate(4)
      write( nowdate(15:16), '(i2)' ) time_nowdate(5)
      write( nowdate(18:19), '(i2)' ) time_nowdate(6)
      do i = 1, 19
        if( nowdate(i:i) == ' ' ) nowdate(i:i) = '0'
      end do
 
      log_info("DA_param_estimation_update",*) 'result [time, name, rank-var, mean, std]: ', &
        nowdate,                                             &
        trim(estimation_target(v)),                          &
        param_array(1+ensemble_myrank,v),                    &
        average( ensemble_nprocs, param_array(:,v), undef ), &
        stddev( ensemble_nprocs, param_array(:,v), undef )
    end do
 
    return
  end subroutine estimation_log
 
  !---------------------------------------------------------------------
  !  Based on the NICAM-LETKF parameter estimation code
  !  created by by Shunji Kotsuki, RIKEN AICS (Dec 2015)
  !---------------------------------------------------------------------
  function phys2func(pidx,phys) ! phys --> func
    implicit none
  
    real(rp) :: phys2func
 
    integer,  intent(in) :: pidx
    real(rp), intent(in) :: phys
 
    real(rp) :: pmax, pmin, pmean
    real(rp) :: tmp_phys
  
    pmax = pest_ulimit(pidx)
    pmin = pest_llimit(pidx)
  
    pmean = 0.5_rp * (pmax + pmin)
    tmp_phys = (phys - pmean) / (pmax - pmean)
  
    phys2func = 0.5_rp * log( (1.0_rp + tmp_phys) / (1.0_rp - tmp_phys) )
    !phys2func = tanh((phys - pmean) / (pmax - pmean))
 
    if (phys <= pmin .or. pmax <= phys) then
      log_info('DA_param_estimation (phys2func)',*) 'Waring: Estimated parameter value is outside of the limit! (pmin/pmax/phys)', pmin, pmax, phys
    endif
  
    return
  end function phys2func
 
  function func2phys(pidx,func) ! func --> phys
    implicit none
  
    real(rp) :: func2phys
 
    integer,  intent(in) :: pidx
    real(rp), intent(in) :: func
 
    real(rp) :: pmax, pmin, pmean
  
    pmax = pest_ulimit(pidx)
    pmin = pest_llimit(pidx)
  
    pmean = 0.5_rp * (pmax + pmin)
  
    func2phys = tanh(func) * (pmax - pmean) + pmean
    !func2phys = 0.5_RP * log( (1.0_RP + func) / (1.0_RP - func) ) * (pmax - pmean) + pmean
  
    if (func2phys <= pmin .or. pmax <= func2phys) then
      log_info('DA_param_estimation (func2phys)',*) 'Waring: Estimated parameter value is outside of the limit! (pmin/pmax/phys)', pmin, pmax, func2phys
    endif
 
    return
  end function func2phys
 
end module mod_da_param_estimation