anchor.diff:    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-14            ! minimum value
  -  real(8), private, parameter :: eps_r = 1.D-9           ! minimum value (reduce version)
  +  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
  
positive.code:     use mod_misc, only: &
         MISC_get_latlon
      use mod_adm, only: &
         ADM_KNONE
      use mod_grd, only: &
         GRD_x,          &
         GRD_XDIR,       &
         GRD_YDIR,       &
         GRD_ZDIR,       &
         GRD_vz
      use mod_cnst, only :  &
         CNST_RAIR, &
           CNST_EGRAV
      use mod_runconf, only: &
         TRC_vmax
      implicit none
  
negative.code:     integer, intent(in)    :: lall
      real(8), intent(inout) :: DIAG_var(ijdim,kdim,lall,6+TRC_VMAX)
  
      real(8) :: pre(kdim), tem(kdim), dz(kdim)
      real(8) :: pre_sfc, tem_sfc
      real(8) :: pre_save
  
      real(8) :: dT, f, df
      real(8) :: lat, lon
  
      integer :: n, k, l, itr
      !---------------------------------------------------------------------------
  
      DIAG_var(:,:,:,:) = 0.D0
  
      do l = 1, lall
      do n = 1, ijdim
  
         dz(ADM_kmin) = GRD_vz(n,ADM_kmin,l,GRD_Z)
         do k = ADM_kmin+1, ADM_kmax+1
            dz(k) = GRD_vz(n,k,l,GRD_Z) - GRD_vz(n,k-1,l,GRD_Z)
         enddo
  
         call MISC_get_latlon( lat, lon,               &
                               GRD_x(n,ADM_KNONE,l,GRD_XDIR), &
                               GRD_x(n,ADM_KNONE,l,GRD_YDIR), &
                               GRD_x(n,ADM_KNONE,l,GRD_ZDIR)  )
  
         pre_sfc = PRE00
  !       tem_sfc = 300.D0
         tem_sfc = 315.D0 - deltaT*sin(lat)**2
  
         !---< from ground surface to lowermost atmosphere >---
         k = ADM_kmin
         ! first guess
         pre(k) = pre_sfc
         tem(k) = tem_sfc
  
         ! Newton-Lapson
         do itr = 1, itrmax
            pre_save = pre(k) ! save
  
            f  = log(pre(k)/pre_sfc) / dz(k) + GRAV / ( RAIR * 0.5D0 * (tem(k)+tem_sfc) )
            df = 1.D0 / (pre(k)*dz(k))
  
            pre(k) = pre(k) - f / df
  !          tem(k) = 300.D0 * ( pre(k)/PRE00 )**KAPPA
            tem(k) = ( 315.D0 - deltaT*sin(lat)**2 - deltaTh*log(pre(k)/PRE00)*cos(lat)**2 ) &
                   * ( pre(k)/PRE00 )**KAPPA
            tem(k) = max( 200.D0, tem(k) )
  
            if( abs(pre_save-pre(k)) <= eps ) exit
         enddo
  
         if ( itr > itrmax ) then
            write(ADM_LOG_FID,*) 'xxx iteration not converged!', k, pre_save-pre(k), pre(k), pre_sfc, tem(k), tem_sfc
            write(*,          *) 'xxx iteration not converged!', k, pre_save-pre(k), pre(k), pre_sfc, tem(k), tem_sfc
            stop
         endif
  
         !---< from lowermost to uppermost atmosphere >---
         do k = ADM_kmin+1, ADM_kmax+1
  
            ! first guess
            pre(k) = pre(k-1)
            tem(k) = 300.D0 * ( pre(k)/PRE00 )**KAPPA
            tem(k) = max( 200.D0, tem(k) )
  
            ! Newton-Lapson
            do itr = 1, itrmax
               pre_save = pre(k) ! save
  
               f  = log(pre(k)/pre(k-1)) / dz(k) + GRAV / ( RAIR * 0.5D0 * (tem(k)+tem(k-1)) )
               df = 1.D0 / (pre(k)*dz(k))
  
               pre(k) = pre(k) - f / df
  !             tem(k) = 300.D0 * ( pre(k)/PRE00 )**KAPPA
               tem(k) = ( 315.D0 - deltaT*sin(lat)**2 - deltaTh*log(pre(k)/PRE00)*cos(lat)**2 ) &
                      * ( pre(k)/PRE00 )**KAPPA
               tem(k) = max( 200.D0, tem(k) )
  
               if( abs(pre_save-pre(k)) <= eps ) exit
  
            enddo
  
            if ( itr > itrmax ) then
               write(ADM_LOG_FID,*) 'xxx iteration not converged!', k, pre_save-pre(k), pre(k), pre(k-1), tem(k), tem(k-1)
               write(*,          *) 'xxx iteration not converged!', k, pre_save-pre(k), pre(k), pre(k-1), tem(k), tem(k-1)
               stop
            endif
         enddo
  
         DIAG_var(n,ADM_kmin-1,l,1) = pre_sfc ! tentative
         DIAG_var(n,ADM_kmin-1,l,2) = tem_sfc ! tentative
         do k = ADM_kmin, ADM_kmax+1
            DIAG_var(n,k,l,1) = pre(k)
            DIAG_var(n,k,l,2) = tem(k)
         enddo
  
      enddo
      enddo
  
      return
    end subroutine hs_init  
  
  