anchor.diff:      use mod_misc, only: &
          MISC_get_latlon
       use mod_adm, only: &
  -       ADM_KNONE
  +       ADM_KNONE, &
  +       ADM_kmin,  &
  +       ADM_kmax
  +    use mod_cnst, only :  &
  +       GRAV  => CNST_EGRAV, &
  +       RAIR  => CNST_RAIR,  &
  +       KAPPA => CNST_KAPPA,  &
  +       PRE00 => CNST_PRE00
       use mod_grd, only: &
  -       GRD_x,          &
  -       GRD_XDIR,       &
  -       GRD_YDIR,       &
  -       GRD_ZDIR,       &
  -       GRD_vz
  -    use mod_cnst, only :  &
  -       CNST_RAIR, &
  -         CNST_EGRAV
  +       GRD_x,    &
  +       GRD_XDIR, &
  +       GRD_YDIR, &
  +       GRD_ZDIR, &
  +       GRD_vz,   &
  +       GRD_Z
       use mod_runconf, only: &
          TRC_vmax
       implicit none
  
positive.code:     integer, intent(in)    :: lall
      real(8), intent(inout) :: DIAG_var(ijdim,kdim,lall,6+TRC_VMAX)
  
      real(8), parameter :: p_surf = 1.D+5         ! [Pa]
      real(8), parameter :: LASPdry = 9.771958146487295D-3
      real(8) :: t_1stlev, t_surf, t_atthelev, t_guess
      real(8) :: p_1stlev(2)
      real(8) :: dZ, dT, f, df
      real(8) :: lat, lon, e
  
      integer :: n, k, l, itr, K0
      !---------------------------------------------------------------------------
  
      K0 = ADM_KNONE
      e = exp(1.D0)
  
      DIAG_var(:,:,:,:) = 0.D0
  
      do l = 1, lall
      do n = 1, ijdim
         k = 1
         call MISC_get_latlon( lat, lon,               &
                               GRD_x(n,K0,l,GRD_XDIR), &
                               GRD_x(n,K0,l,GRD_YDIR), &
                               GRD_x(n,K0,l,GRD_ZDIR)  )
         t_surf = 315.D0 - deltaT*sin(lat)**2.D0
         p_1stlev(:) = p_surf
         dZ = GRD_vz(n,k,l,1) - 0.D0
  
         ! Newton-Lapson
         p_1stlev(1) = 0.D0 ! to do the first step
         do itr = 1, itrmax
            if( abs(p_1stlev(1)-p_1stlev(2)) <= eps_r ) exit
  
            p_1stlev(1) = p_1stlev(2)
            t_1stlev = ( 315.D0 - deltaT*sin(lat)**2 - deltaTh*log(p_1stlev(1)/p0) &
                                              * cos(lat)**2 ) * (p_1stlev(1)/p0)**kai
            f  = ( log( p_surf ) - log( p_1stlev(1) ) ) / ( dZ ) &
               + CNST_EGRAV / CNST_RAIR * 2.D0 / ( t_surf + t_1stlev )
            df = 1.D0 / dZ / p_1stlev(1)
  
            p_1stlev(2) = p_1stlev(1) + f/df
         enddo
         if ( itr > itrmax ) then
            write(*,*) 'xxx iteration not converged!', &
                       p_1stlev(1)-p_1stlev(2), p_1stlev(2), t_1stlev, t_surf
         endif
  
         DIAG_var(n,k,l,1) = p_1stlev(2)
         t_atthelev = ( 315.D0 - deltaT*sin(lat)**2 - deltaTh*log(DIAG_var(n,k,l,1)/p0) &
                                           * cos(lat)**2 ) * (DIAG_var(n,k,l,1)/p0)**kai
         DIAG_var(n,k,l,2) = max(200.D0, t_atthelev)
  
         ! guess pressure field upper k=1
         do k=2, kdim
            dZ = GRD_vz(n,k,l,1) - GRD_vz(n,k-1,l,1)
            t_guess = DIAG_var(n,k-1,l,2) - LASPdry * dZ
            dT = ( t_guess + DIAG_var(n,k-1,l,2) )/2.D0
            DIAG_var(n,k,l,1) = DIAG_var(n,k-1,l,1) * e**( -1.D0 * dZ * (g/(R*dT)) )
            t_atthelev = ( 315.D0 - deltaT*sin(lat)**2 - deltaTh*log(DIAG_var(n,k,l,1)/p0) &
                                                     * cos(lat)**2 ) * (DIAG_var(n,k,l,1)/p0)**kai
            DIAG_var(n,k,l,2) = max(200.D0, t_atthelev)
         enddo
  
      enddo
      enddo
   
      return
    end subroutine hs_init  
  
  
negative.code:   real(8), private, parameter :: pi = 3.14159265358979323846D0
    real(8), private, parameter :: d2r = pi/180.D0          ! Degree to Radian
    real(8), private, parameter :: r2d = 180.D0/pi          ! Radian to Degree
    real(8), private, parameter :: eps = 1.D-8              ! minimum value
    real(8), private, parameter :: zero = 0.D0              ! zero
  
    ! test configurations
    integer, private, parameter :: PRCS_D = 8
    ! for Held and Suarez
  