anchor.diff:    end subroutine bsstate_setup
   
     !-----------------------------------------------------------------------------
  -  subroutine output_info
  -    use mod_adm, only: &
  -       ADM_LOG_FID, &
  -       ADM_kall,    &
  -       ADM_kmin,    &
  -       ADM_kmax
  -    implicit none
  -
  -    integer :: k
  -    !---------------------------------------------------------------------------
  -
  -    write(ADM_LOG_FID,*) 'Msg : Sub[nhm_bs_output_info]/Mod[basicstate]'
  -    write(ADM_LOG_FID,*) ' --- Basic state type         : ', trim(ref_type)
  -    write(ADM_LOG_FID,*) ' --- Reference pressure       : ', pre_g
  -    write(ADM_LOG_FID,*) ' --- Reference temperature    : ', tem_g
  -    write(ADM_LOG_FID,*) ' --- Reference density        : ', rho_g
  -    write(ADM_LOG_FID,*) ' --- Vaisala frequency        : ', BV_freq
  -    write(ADM_LOG_FID,*) ' --- Lapse rate of temperature: ', TGAMMA
  -    write(ADM_LOG_FID,*) ' --- Effective height         : ', ZT
  -
  -    write(ADM_LOG_FID,*) '-------------------------------------------------------'
  -    write(ADM_LOG_FID,*) 'Level   Density  Pressure     Temp. Pot. Tem.        qv'
  -    do k=1,ADM_kall
  -       write(ADM_LOG_FID,'(I4,F12.4,F10.2,F10.2,F10.2,F10.7)') &
  -                         k, rho_ref(k), pre_ref(k), tem_ref(k),th_ref(k), qv_ref(k)
  -       if (  k == ADM_kmin-1 ) write(ADM_LOG_FID,*) '-------------------------------------------------------'
  -       if (  k == ADM_kmax   ) write(ADM_LOG_FID,*) '-------------------------------------------------------'
  -    enddo
  -
  -    return
  -  end subroutine output_info
  -
  -  !-----------------------------------------------------------------------------
  -  subroutine set_referencestate
  -    use mod_cnst, only : &
  -         CNST_EGRAV,     &
  -         CNST_RAIR,      &
  -         CNST_RVAP,      &
  -         CNST_PRE00,     &
  -         CNST_KAPPA
  -    use mod_adm, only :  &
  -       ADM_have_pl, &
  -         ADM_kall,       &
  -         ADM_kmin,       &
  -         ADM_kmax
  -    use mod_grd, only :  &
  -         GRD_gz,         &
  -         GRD_dgz,        &
  -         GRD_gzh,        &
  -         GRD_afac,       &
  -         GRD_bfac
  -    implicit none
  -
  -    real(RP) :: dpre_ref_k
  -    real(RP) :: pre_s, rho_s, total_mass0, total_mass, mass_diff_ratio
  -
  -    integer :: k
  -    !---------------------------------------------------------------------------
  -
  -    !--- calculation of reference geopotential
  -    do k=1,ADM_kall
  -       phi_ref(k) = CNST_EGRAV * GRD_gz(k)
  -    enddo
  -    !
  -    if ( ref_type=='NOBASE') then
  -       !
  -       phi_ref = 0.0_RP
  -       rho_ref = 0.0_RP
  -       pre_ref = 0.0_RP
  -       tem_ref = 0.0_RP
  -       th_ref  = 0.0_RP
  -       qv_ref  = 0.0_RP
  -       !
  -       pre_bs    = 0.0_RP
  -       pre_bs_pl = 0.0_RP
  -       tem_bs    = 0.0_RP
  -       tem_bs_pl = 0.0_RP
  -       th_bs    = 0.0_RP
  -       th_bs_pl = 0.0_RP
  -       qv_bs    = 0.0_RP
  -       qv_bs_pl = 0.0_RP
  -       rho_bs    = 0.0_RP
  -       rho_bs_pl = 0.0_RP
  -       !
  -       ! 04/12/25 M.Satoh add
  -       if ( tem_g /= 0.0_RP ) then
  -          rho_g = pre_g / CNST_RAIR / tem_g
  -       else
  -          rho_g = 0.0_RP
  -       endif
  -       !
  -    else if ( ref_type == 'TEM' ) then
  -       qv_ref = 0.0_RP
  -       !
  -       !---  calculation of reference temperature
  -       do k=1,ADM_kall
  -          if ( GRD_gz(k) <=ZT ) then
  -             tem_ref(k) = tem_g - TGAMMA * GRD_gz(k)
  -          else
  -             tem_ref(k) = tem_g - TGAMMA * ZT
  -          endif
  -       enddo
  -       !
  -       !--- calculation of density at the surface
  -       rho_g = pre_g / CNST_RAIR / tem_g
  -       !
  -       !--- calculation of reference pressure and density
  -       !--- just below the ground level
  -       pre_ref(ADM_kmin-1)                                &
  -            = pre_g                                        &
  -            + 0.5_RP                                       &
  -            * ( phi_ref(ADM_kmin) - phi_ref(ADM_kmin-1) ) &
  -            * rho_g
  -       rho_ref(ADM_kmin-1)        &
  -            = pre_ref(ADM_kmin-1) &
  -            / CNST_RAIR           &
  -            / tem_ref(ADM_kmin-1)
  -       !
  -       !--- Reference pressure and density at the first level
  -       pre_ref(ADM_kmin)                                        &
  -            = pre_ref(ADM_kmin-1)                               &
  -            - ( phi_ref(ADM_kmin) - phi_ref(ADM_kmin-1) ) * rho_g
  -       rho_ref(ADM_kmin)        &
  -            = pre_ref(ADM_kmin) &
  -            / CNST_RAIR         &
  -            / tem_ref(ADM_kmin)
  -       !
  -       !--- Reference pressure and density at k level
  -       !---    In this caluculation, the hydrostatic balance equation
  -       !---    ( dP/dz=-RHO dPHI/dz ) is applied at the half integer
  -       !---    level ( k-1/2 ). RHO is obtained by extrapolation from
  -       !---    the values at lower level.
  -       !---
  -       !--- Consistent way(?) in the scheme.
  -       !---
  -       !--- pre_ref(k) - pre_ref(k-1)
  -       !--- = - 0.5_RP * ( GRD_afac(k) * rho_ref(k)
  -       !---              +GRD_bfac(k) * rho_ref(k-1)  )
  -       !---   * ( phi_ref(k) - phi_ref(k-1) )
  -       !---
  -       !--- rho_ref(k)*CNST_RAIR*tem_ref(k)  - pre_ref(k-1)
  -       !--- = - 0.5_RP * ( GRD_afac(k) * rho_ref(k)
  -       !---              +GRD_bfac(k) * rho_ref(k-1)  )
  -       !---   * ( phi_ref(k) - phi_ref(k-1) )
  -       !
  -       !--- rho_ref(k)*( CNST_RAIR*tem_ref(k)
  -       !---            + 0.5_RP * GRD_afac(k)
  -       !---            * ( phi_ref(k) - phi_ref(k-1) )
  -       !---            )
  -       !--- = pre_ref(k-1)
  -       !---            - 0.5_RP * GRD_bfac(k) * rho_ref(k-1)
  -       !---            * ( phi_ref(k) - phi_ref(k-1) )
  -       !---
  -       do k = ADM_kmin+1, ADM_kmax+1
  -          rho_ref(k) = &
  -               ( pre_ref(k-1) &
  -               - 0.5_RP * GRD_bfac(k) * rho_ref(k-1) &
  -               * ( phi_ref(k) - phi_ref(k-1) )&
  -               ) / &
  -               ( CNST_RAIR*tem_ref(k) &
  -               + 0.5_RP * GRD_afac(k)  &
  -               * ( phi_ref(k) - phi_ref(k-1) )  &
  -               )
  -          pre_ref(k) = rho_ref(k) * CNST_RAIR * tem_ref(k)
  -       enddo
  -       !
  -       !--- calculation of reference potential temperature
  -       do k = 1, ADM_kall
  -          th_ref(k)                                        &
  -               = tem_ref(k)                                &
  -               * ( CNST_PRE00 / pre_ref(k) ) ** CNST_KAPPA
  -       enddo
  -       !
  -    else if ( ref_type == 'RHO' ) then
  -       qv_ref = 0.0_RP
  -       !
  -       !---  calculation of reference density
  -       total_mass0 = CNST_PRE00/CNST_EGRAV
  -       pre_s = CNST_PRE00
  -       do
  -          do k = ADM_kmin-1, ADM_kmax+1
  -             rho_ref(k) = pre_s/CNST_RAIR/tem_g / exp(CNST_EGRAV*GRD_gz(k)/CNST_RAIR/tem_g)
  -          enddo
  -          total_mass = 0.0_RP
  -          do k = ADM_kmin, ADM_kmax
  -             total_mass = total_mass + GRD_dgz(k)*rho_ref(k)
  -          enddo
  -          mass_diff_ratio = total_mass0/total_mass
  -          if ( abs(mass_diff_ratio-1.0_RP)<1.E-8_RP) then
  -             exit
  -          else
  -             pre_s = pre_s * mass_diff_ratio
  -          endif
  -       enddo
  -       !
  -       !--- calculation of density at the surface
  -       rho_s = pre_s /CNST_RAIR / tem_g
  -       !
  -       !--- calculation of reference pressure and density
  -       !--- just below the ground level
  -       pre_ref(ADM_kmin-1)                                &
  -            = pre_s                                       &
  -            + 0.5_RP                                       &
  -            * ( phi_ref(ADM_kmin) - phi_ref(ADM_kmin-1) ) &
  -            * rho_s
  -       !
  -       !--- Reference pressure and density at the first level
  -       pre_ref(ADM_kmin)                                        &
  -            = pre_ref(ADM_kmin-1)                               &
  -            - ( phi_ref(ADM_kmin) - phi_ref(ADM_kmin-1) ) * rho_s
  -       !
  -       !--- Reference pressure and density at k level
  -       do k = ADM_kmin+1, ADM_kmax+1
  -          pre_ref(k) = pre_ref(k-1) &
  -               - 0.5_RP * ( GRD_afac(k) * rho_ref(k)     &
  -                          +GRD_bfac(k) * rho_ref(k-1)  )&
  -               * ( phi_ref(k) - phi_ref(k-1) )
  -       enddo
  -       !
  -       !--- calculation of reference temperature & pot temperature
  -       do k = 1, ADM_kall
  -          tem_ref(k) = pre_ref(k) / rho_ref(k) / CNST_RAIR
  -          th_ref(k)                                        &
  -               = tem_ref(k)                                &
  -               * ( CNST_PRE00 / pre_ref(k) ) ** CNST_KAPPA
  -       enddo
  -       !
  -    else if ( ref_type == 'TH' ) then
  -       qv_ref = 0.0_RP
  -       !
  -       !--- calculation of reference pot. temp.
  -       !--- just below the surface level.
  -       th_ref(ADM_kmin-1)                                 &
  -            = th_g                                         &
  -            / exp ( BV_freq**2 / CNST_EGRAV                     &
  -                   * ( GRD_gzh(ADM_kmin) - GRD_gz(ADM_kmin-1) ) )
  -       !
  -       !--- calculation of reference pot. temp. at the first level.
  -       th_ref(ADM_kmin)                                   &
  -            = th_g                                         &
  -            * exp ( BV_freq**2 / CNST_EGRAV                     &
  -                   * ( GRD_gz(ADM_kmin) - GRD_gzh(ADM_kmin) ) )
  -       !
  -       !--- calculation of reference pot. temp. at k level
  -       do k = ADM_kmin+1, ADM_kmax+1
  -          if ( GRD_gz(k) <= ZT ) then
  -             th_ref(k) = th_ref(k-1)                     &
  -                  * exp ( BV_freq**2 / CNST_EGRAV              &
  -                         * ( GRD_gz(k) - GRD_gz(k-1) ) )
  -          else
  -             th_ref(k) = th_ref(k-1)
  -          endif
  -       enddo
  -       !
  -       !--- calculation of density at the surface
  -       tem_g = th_g / ( CNST_PRE00 / pre_g ) ** CNST_KAPPA
  -       rho_g = pre_g / CNST_RAIR / tem_g
  -       !
  -       !--- calculation of reference pressure, temperature, density
  -       !--- just below the ground level.
  -       pre_ref(ADM_kmin-1)                                &
  -            = pre_g                                        &
  -            + 0.5_RP                                       &
  -            * ( phi_ref(ADM_kmin) - phi_ref(ADM_kmin-1) ) &
  -            * rho_g
  -       tem_ref(ADM_kmin-1)                                        &
  -            = th_ref(ADM_kmin-1)                                  &
  -            / ( CNST_PRE00 / pre_ref(ADM_kmin-1) ) ** CNST_KAPPA
  -       rho_ref(ADM_kmin-1)        &
  -            = pre_ref(ADM_kmin-1) &
  -            / CNST_RAIR           &
  -            / tem_ref(ADM_kmin-1)
  -       !
  -       !--- calculation of reference pressure and density
  -       !--- at the first level
  -       pre_ref(ADM_kmin)                                         &
  -            = pre_ref(ADM_kmin-1)                                &
  -            - ( phi_ref(ADM_kmin) - phi_ref(ADM_kmin-1) ) * rho_g
  -       tem_ref(ADM_kmin)                                       &
  -            = th_ref(ADM_kmin)                                 &
  -            / ( CNST_PRE00 / pre_ref(ADM_kmin) ) ** CNST_KAPPA
  -       rho_ref(ADM_kmin)                                       &
  -            = pre_ref(ADM_kmin) / CNST_RAIR / tem_ref(ADM_kmin)
  -       !
  -       !--- Reference pressure and density at k level
  -       !---    In this caluculation, the hydrostatic balance equation
  -       !---    ( dP/dz=-RHO dPHI/dz ) is applied at the half integer
  -       !---    level ( k-1/2 ). RHO is obtained by extrapolation from
  -       !---    the values at lower level.
  -       !
  -       !--- fist guess
  -       do k = ADM_kmin+1, ADM_kmax+1
  -          pre_ref(k) = pre_ref(k-1)              &
  -               - ( phi_ref(k) - phi_ref(k-1) )   &
  -               * ( rho_ref(k-1)                  &
  -               + ( rho_ref(k-1) - rho_ref(k-2) ) &
  -               / ( GRD_gz(k-1)  - GRD_gz(k-2) )   &
  -               * ( GRD_gz(k)    - GRD_gz(k-1) ) * 0.5_RP )
  -          tem_ref(k)                                      &
  -               = th_ref(k)                                &
  -               / ( CNST_PRE00 / pre_ref(k) ) ** CNST_KAPPA
  -          rho_ref(k)        &
  -               = pre_ref(k) &
  -               / CNST_RAIR  &
  -               / tem_ref(k)
  -       enddo
  -       !--- hydro static balance adjustment
  -       do k = ADM_kmin+1, ADM_kmax+1
  -          do
  -             tem_ref(k) = th_ref(k)                                &
  -                  / ( CNST_PRE00 / pre_ref(k) ) ** CNST_KAPPA
  -             rho_ref(k) = &
  -                  ( pre_ref(k-1) &
  -                  - 0.5_RP * GRD_bfac(k) * rho_ref(k-1) &
  -                  * ( phi_ref(k) - phi_ref(k-1) )&
  -                  ) / &
  -                  ( CNST_RAIR*tem_ref(k) &
  -                  + 0.5_RP * GRD_afac(k)  &
  -                  * ( phi_ref(k) - phi_ref(k-1) )  &
  -                  )
  -             dpre_ref_k = rho_ref(k) * CNST_RAIR * tem_ref(k) - pre_ref(k)
  -             pre_ref(k) = pre_ref(k)+ dpre_ref_k
  -             if ( abs(dpre_ref_k) < 1.E-10_RP ) exit
  -             if ( ADM_have_pl ) then
  -                write(*,*) k,abs(dpre_ref_k)
  -             endif
  -          enddo
  -       enddo
  -    else if ( ref_type == 'TH-SP' ) then
  -       qv_ref = 0.0_RP
  -       !
  -       !--- calculation of reference pot. temp.
  -       !--- just below the surface level.
  -       do k=ADM_kmin-1,ADM_kmax+1
  -          if ( GRD_gz(k)<10000.0_RP) then
  -             th_ref(k) = th_g
  -          else
  -             th_ref(k) = th_g+10.0_RP/1000.0_RP*(GRD_gz(k)-10000.0_RP)
  -          endif
  -       enddo
  -       !
  -       !--- calculation of density at the surface
  -       tem_g = th_g / ( CNST_PRE00 / pre_g ) ** CNST_KAPPA
  -       rho_g = pre_g / CNST_RAIR / tem_g
  -       !
  -       !--- calculation of reference pressure, temperature, density
  -       !--- just below the ground level.
  -       pre_ref(ADM_kmin-1)                                &
  -            = pre_g                                        &
  -            + 0.5_RP                                       &
  -            * ( phi_ref(ADM_kmin) - phi_ref(ADM_kmin-1) ) &
  -            * rho_g
  -       tem_ref(ADM_kmin-1)                                        &
  -            = th_ref(ADM_kmin-1)                                  &
  -            / ( CNST_PRE00 / pre_ref(ADM_kmin-1) ) ** CNST_KAPPA
  -       rho_ref(ADM_kmin-1)        &
  -            = pre_ref(ADM_kmin-1) &
  -            / CNST_RAIR           &
  -            / tem_ref(ADM_kmin-1)
  -       !
  -       !--- calculation of reference pressure and density
  -       !--- at the first level
  -       pre_ref(ADM_kmin)                                         &
  -            = pre_ref(ADM_kmin-1)                                &
  -            - ( phi_ref(ADM_kmin) - phi_ref(ADM_kmin-1) ) * rho_g
  -       tem_ref(ADM_kmin)                                       &
  -            = th_ref(ADM_kmin)                                 &
  -            / ( CNST_PRE00 / pre_ref(ADM_kmin) ) ** CNST_KAPPA
  -       rho_ref(ADM_kmin)                                       &
  -            = pre_ref(ADM_kmin) / CNST_RAIR / tem_ref(ADM_kmin)
  -       !
  -       !--- Reference pressure and density at k level
  -       !---    In this caluculation, the hydrostatic balance equation
  -       !---    ( dP/dz=-RHO dPHI/dz ) is applied at the half integer
  -       !---    level ( k-1/2 ). RHO is obtained by extrapolation from
  -       !---    the values at lower level.
  -       !
  -       !--- fist guess
  -       do k = ADM_kmin+1, ADM_kmax+1
  -          pre_ref(k) = pre_ref(k-1)              &
  -               - ( phi_ref(k) - phi_ref(k-1) )   &
  -               * ( rho_ref(k-1)                  &
  -               + ( rho_ref(k-1) - rho_ref(k-2) ) &
  -               / ( GRD_gz(k-1)  - GRD_gz(k-2) )   &
  -               * ( GRD_gz(k)    - GRD_gz(k-1) ) * 0.5_RP )
  -          tem_ref(k)                                      &
  -               = th_ref(k)                                &
  -               / ( CNST_PRE00 / pre_ref(k) ) ** CNST_KAPPA
  -          rho_ref(k)        &
  -               = pre_ref(k) &
  -               / CNST_RAIR  &
  -               / tem_ref(k)
  -       enddo
  -       !--- hydro static balance adjustment
  -       do k = ADM_kmin+1, ADM_kmax+1
  -          do
  -             tem_ref(k) = th_ref(k)                                &
  -                  / ( CNST_PRE00 / pre_ref(k) ) ** CNST_KAPPA
  -             rho_ref(k) = &
  -                  ( pre_ref(k-1) &
  -                  - 0.5_RP * GRD_bfac(k) * rho_ref(k-1) &
  -                  * ( phi_ref(k) - phi_ref(k-1) )&
  -                  ) / &
  -                  ( CNST_RAIR*tem_ref(k) &
  -                  + 0.5_RP * GRD_afac(k)  &
  -                  * ( phi_ref(k) - phi_ref(k-1) )  &
  -                  )
  -             dpre_ref_k = rho_ref(k) * CNST_RAIR * tem_ref(k) - pre_ref(k)
  -             pre_ref(k) = pre_ref(k)+ dpre_ref_k
  -             if ( abs(dpre_ref_k) < 1.E-10_RP ) exit
  -             if ( ADM_have_pl ) then
  -                write(*,*) k,abs(dpre_ref_k)
  -             endif
  -          enddo
  -       enddo
  -    else if ( ref_type == 'OOYAMA' ) then
  -       call ooyama_reference
  -       do k = 1, ADM_kall
  -          pre_ref(k) = rho_ref(k) * tem_ref(k) &
  -               * ( (1.0_RP-qv_ref(k))*CNST_RAIR+qv_ref(k)*CNST_RVAP )
  -          th_ref(k)                                        &
  -               = tem_ref(k)                                &
  -               * ( CNST_PRE00 / pre_ref(k) ) ** CNST_KAPPA
  -       enddo
  -    endif
  -
  -!   if ( ref_type == 'GCSS_CASE1' ) then
  -    if ( ref_type == 'GCSS_CASE1' .or. ref_type=='TWP-ICE') then ! 11/08/13 A.Noda [add]
  -       call gcss_reference
  -       do k = 1, ADM_kall
  -          tem_ref(k) = th_ref(k)                                &
  -               / ( CNST_PRE00 / pre_ref(k) ) ** CNST_KAPPA
  -          rho_ref(k) = pre_ref(k)/tem_ref(k) &
  -               / ( (1.0_RP-qv_ref(k))*CNST_RAIR+qv_ref(k)*CNST_RVAP )
  -       enddo
  -    endif
  -
  -
  -  end subroutine set_referencestate
  -
  -  !-----------------------------------------------------------------------------
  -  subroutine gcss_reference
  -    use mod_misc, only: &
  -       MISC_get_available_fid
  -    implicit none
  -
  -    integer :: fid
  -
  -    fid = MISC_get_available_fid()
  -    Open(fid,file=Trim(ref_fname),status='old',form='unformatted')
  -    read(fid) th_ref(:)
  -    read(fid) pre_ref(:)
  -    read(fid) qv_ref(:)
  -    close(fid)
  -
  -  end subroutine gcss_reference
  -
  -  !-----------------------------------------------------------------------------
  -  subroutine ooyama_reference
  +  subroutine bsstate_input_ref( basename )
       use mod_misc, only: &
          MISC_get_available_fid
       use mod_adm, only: &
  -       ADM_kmin, &
  -       ADM_kmax
  +       ADM_kall
  +    use mod_cnst, only: &
  +       GRAV  => CNST_EGRAV, &
  +       Rdry  => CNST_RAIR,  &
  +       CPdry => CNST_CP,    &
  +       CVdry => CNST_CV,    &
  +       Rvap  => CNST_RVAP,  &
  +       PRE00 => CNST_PRE00
       use mod_grd, only: &
          GRD_gz
       implicit none
   
  -    character(len=128) :: fname = 'ooyama_profile.dat'
  -
  -    real(RP), allocatable :: z_s  (:)
  -    real(RP), allocatable :: rho_s(:)
  -    real(RP), allocatable :: tem_s(:)
  -    real(RP), allocatable :: qv_s (:)
  +    Character(*), Intent(in) :: basename
  +
  +    real(RP) :: kappa
   
       integer :: fid
  -    integer :: kmax
  -    integer :: k, kk, kp
  -
  -    real(RP) :: lag_intpl
  -    real(RP) :: z,z1,p1,z2,p2,z3,p3
  -    lag_intpl(z,z1,p1,z2,p2,z3,p3) = ((z-z2)*(z-z3))/((z1-z2)*(z1-z3))*p1 &
  -                                   + ((z-z1)*(z-z3))/((z2-z1)*(z2-z3))*p2 &
  -                                   + ((z-z1)*(z-z2))/((z3-z1)*(z3-z2))*p3
  +    integer :: k
       !---------------------------------------------------------------------------
   
  -    !--- read sounding data ( ooyama(2001) )
  +    kappa = Rdry / CPdry
  +
       fid = MISC_get_available_fid()
  -    open( unit   = fid,         &
  -          file   = trim(fname), &
  -          status = 'old',       &
  -          form   = 'unformatted')
  -
  -       read(fid) kmax
  -
  -       allocate( z_s  (kmax) )
  -       allocate( rho_s(kmax) )
  -       allocate( tem_s(kmax) )
  -       allocate( qv_s (kmax) )
  -
  -       read(fid) z_s  (:)
  -       read(fid) rho_s(:)
  -       read(fid) tem_s(:)
  -       read(fid) qv_s (:)
  +    open( unit   = fid,            &
  +          file   = trim(basename), &
  +          status = 'old',          &
  +          form   = 'unformatted'   )
  +       read(fid) pre_ref(:)
  +       read(fid) tem_ref(:)
  +       read(fid) qv_ref (:)
       close(fid)
   
  -    !--- initialization
  -    do k = ADM_kmin, ADM_kmax
  -    do kk = 1, kmax-1
  -       if (       z_s(kk)   <= GRD_gz(k) &
  -            .AND. z_s(kk+1) >= GRD_gz(k) ) then
  -          if ( kk == 1 ) then
  -             kp = 2
  -          else
  -             kp = kk
  -          endif
  -
  -          rho_ref(k) = lag_intpl(GRD_gz(k),z_s(kp+1),rho_s(kp+1),z_s(kp),rho_s(kp),z_s(kp-1),rho_s(kp-1))
  -          tem_ref(k) = lag_intpl(GRD_gz(k),z_s(kp+1),tem_s(kp+1),z_s(kp),tem_s(kp),z_s(kp-1),tem_s(kp-1))
  -          qv_ref (k) = lag_intpl(GRD_gz(k),z_s(kp+1),qv_s (kp+1),z_s(kp),qv_s (kp),z_s(kp-1),qv_s (kp-1))
  -          exit
  -       endif
  -    enddo
  -    enddo
  -
  -    k  = ADM_kmin-1
  -    kp = 2
  -
  -    rho_ref(k) = lag_intpl(GRD_gz(k),z_s(kp+1),rho_s(kp+1),z_s(kp),rho_s(kp),z_s(kp-1),rho_s(kp-1))
  -    tem_ref(k) = lag_intpl(GRD_gz(k),z_s(kp+1),tem_s(kp+1),z_s(kp),tem_s(kp),z_s(kp-1),tem_s(kp-1))
  -    qv_ref (k) = lag_intpl(GRD_gz(k),z_s(kp+1),qv_s (kp+1),z_s(kp),qv_s (kp),z_s(kp-1),qv_s (kp-1))
  -
  -    k  = ADM_kmax+1
  -    kp = kmax-1
  -
  -    rho_ref(k) = lag_intpl(GRD_gz(k),z_s(kp+1),rho_s(kp+1),z_s(kp),rho_s(kp),z_s(kp-1),rho_s(kp-1))
  -    tem_ref(k) = lag_intpl(GRD_gz(k),z_s(kp+1),tem_s(kp+1),z_s(kp),tem_s(kp),z_s(kp-1),tem_s(kp-1))
  -    qv_ref (k) = lag_intpl(GRD_gz(k),z_s(kp+1),qv_s (kp+1),z_s(kp),qv_s (kp),z_s(kp-1),qv_s (kp-1))
  -
  -    deallocate( rho_s )
  -    deallocate( tem_s )
  -    deallocate( qv_s  )
  +    !--- additional reference state.
  +    do k = 1, ADM_kall
  +       th_ref (k) = tem_ref(k) * ( PRE00 / pre_ref(k) )**kappa
  +       phi_ref(k) = GRAV * GRD_gz(k)
  +       rho_ref(k) = pre_ref(k) / tem_ref(k) / ( ( 1.0_RP - qv_ref(k) ) * Rdry &
  +                                              + (          qv_ref(k) ) * Rvap )
  +    enddo
   
       return
  -  end subroutine ooyama_reference
  +  end subroutine bsstate_input_ref
   
     !-----------------------------------------------------------------------------
     !> generation of basic state from reference state
     subroutine set_basicstate
       use mod_cnst, only: &
  -       CNST_EGRAV
  +       GRAV  => CNST_EGRAV
       use mod_adm, only: &
  -       ADM_have_pl, &
  +       ADM_prc_pl,  &
  +       ADM_prc_me,  &
          ADM_lall,    &
          ADM_lall_pl, &
          ADM_gall,    &
  
positive.code: !-------------------------------------------------------------------------------
  !
  !+  Basic state module
  !
  !-------------------------------------------------------------------------------
  module mod_bsstate
    !-----------------------------------------------------------------------------
    !
    !++ Description:
    !       This module is for the set of basic state for non-hydrostatic
    !       model.
    !
    !
    !++ Current Corresponding Author : H.Tomita
    !
    !++ History:
    !      Version   Date       Comment
    !      -----------------------------------------------------------------------
    !      0.00      04-02-17   Imported from igdc-4.34
    !                11-08-13   A.Noda : add twp-ice exp
    !      -----------------------------------------------------------------------
    !
    !-----------------------------------------------------------------------------
    !
    !++ Used modules
  
negative.code:     write(ADM_LOG_FID,nml=BSSTATEPARAM)
  
      !--- allocation of reference variables
      allocate( phi_ref(ADM_kall) )
      allocate( rho_ref(ADM_kall) )
      allocate( pre_ref(ADM_kall) )
      allocate( tem_ref(ADM_kall) )
      allocate( th_ref (ADM_kall) )
      allocate( qv_ref (ADM_kall) )
      phi_ref(:) = 0.0_DP
      rho_ref(:) = 0.0_DP
      pre_ref(:) = 0.0_DP
      tem_ref(:) = 0.0_DP
      th_ref (:) = 0.0_DP
      qv_ref (:) = 0.0_DP
  
      allocate( rho_bs   (ADM_gall   ,ADM_kall,ADM_lall   ) )
      allocate( rho_bs_pl(ADM_gall_pl,ADM_kall,ADM_lall_pl) )
      rho_bs   (:,:,:) = 0.0_RP
      rho_bs_pl(:,:,:) = 0.0_RP
  
      allocate( pre_bs   (ADM_gall   ,ADM_kall,ADM_lall   ) )
      allocate( pre_bs_pl(ADM_gall_pl,ADM_kall,ADM_lall_pl) )
      pre_bs   (:,:,:) = 0.0_RP
      pre_bs_pl(:,:,:) = 0.0_RP
  
      allocate( tem_bs   (ADM_gall   ,ADM_kall,ADM_lall   ) )
      allocate( tem_bs_pl(ADM_gall_pl,ADM_kall,ADM_lall_pl) )
      tem_bs   (:,:,:) = 0.0_RP
      tem_bs_pl(:,:,:) = 0.0_RP
  
      if ( ref_type == 'INPUT' ) then
         call bsstate_input_ref(ref_fname)
      endif
  
      !--- calculation of basic state
      call set_basicstate
  
      !--- output the information
      call output_info
  
  