anchor.diff:      implicit none
   
       namelist / BSSTATEPARAM / &
  -         ref_type, & !--- type of basic state
  -         ZT,       & !--- if z>ZT, equi-temperature
  -         pre_g,    & !--- reference pressure
  -         tem_g,    & !--- reference temperature
  -         TGAMMA,   & !--- lapse rate when ref_type='TEM'.
  -         th_g,     & !--- reference potential temp. for when ref_type='TH'
  -         BV_freq,  & !--- Vaisala freq when ref_type='TH'
  -         ref_fname   !--- in/output base name if necessary.
  +       ref_type,  &
  +       ref_fname
   
       integer :: ierr
       !---------------------------------------------------------------------------
  
positive.code:   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
      use mod_misc, only: &
         MISC_get_available_fid
      use mod_adm, only: &
         ADM_kmin, &
         ADM_kmax
      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 (:)
  
      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
      !---------------------------------------------------------------------------
  
      !--- read sounding data ( ooyama(2001) )
      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 (:)
      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  )
  
      return
    end subroutine ooyama_reference
  
    !-----------------------------------------------------------------------------
    !> generation of basic state from reference state
    subroutine set_basicstate
      use mod_cnst, only: &
         CNST_EGRAV
      use mod_adm, only: &
         ADM_have_pl, &
         ADM_lall,    &
         ADM_lall_pl, &
         ADM_gall,    &
  
negative.code:   private
    !-----------------------------------------------------------------------------
    !
    !++ Public procedure
    !
    public :: bsstate_setup
  
    !-----------------------------------------------------------------------------
    !
    !++ Public parameters & variables
    !
    real(RP), public, allocatable :: rho_bs   (:,:,:)
    real(RP), public, allocatable :: rho_bs_pl(:,:,:)
    real(RP), public, allocatable :: pre_bs   (:,:,:)
    real(RP), public, allocatable :: pre_bs_pl(:,:,:)
    real(RP), public, allocatable :: tem_bs   (:,:,:)
    real(RP), public, allocatable :: tem_bs_pl(:,:,:)
  
    !-----------------------------------------------------------------------------
    !
    !++ Private procedures
    !
    private :: bsstate_input_ref
    private :: set_basicstate
    private :: output_info
  
    !-----------------------------------------------------------------------------
    !
    !++ Private parameters & variables
    !
    character(len=ADM_NSYS), private :: ref_type = 'NOBASE'  !--- Basic state type
    !                                            = 'NOBASE' : no basic state
    !                                            = 'INPUT'  : input
    !                                            = 'TEM'    : temperature is given.
    !                                            = 'TH'     : potential temperature is given.
  
    character(len=ADM_MAXFNAME), private :: ref_fname = 'ref.dat'
  
    real(DP), private, allocatable :: phi_ref(:) ! reference phi
    real(DP), private, allocatable :: rho_ref(:) ! reference density
    real(DP), private, allocatable :: pre_ref(:) ! reference pressure
    real(DP), private, allocatable :: tem_ref(:) ! reference temperature
    real(DP), private, allocatable :: qv_ref (:) ! water vapor
    real(DP), private, allocatable :: th_ref (:) ! reference potential temperature
  
    !-----------------------------------------------------------------------------
  contains
    !-----------------------------------------------------------------------------
  