anchor.diff:         lat = GMTR_lat(n,l)
          lon = GMTR_lon(n,l)
   
  -!       z_local(ADM_kmin-1) = GRD_vz(n,2,l,GRD_ZH)
  +       z(ADM_kmin-1) = GRD_vz(n,2,l,GRD_ZH)
          do k = ADM_kmin, ADM_kmax+1
  -          z_local = GRD_vz(n,k,l,GRD_Z)
  -
  -          call baroclinic_wave_test( deep,    &  ! [IN ]
  -                                     moist,   &  ! [IN ]
  -                                     pertt,   &  ! [IN ]
  -                                     Xfact,   &  ! [IN ]
  -                                     lon,     &  ! [IN ]
  -                                     lat,     &  ! [IN ]
  -                                     p,       &  ! [IN ]
  -                                     z_local, &  ! [IN ] (z)
  -                                     zcoords, &  ! [IN ]
  +          z(k) = GRD_vz(n,k,l,GRD_Z)
  +       enddo
  +
  +       do k = 1, kdim
  +          call baroclinic_wave_test( deep,       &  ! [IN ]
  +                                     moist,      &  ! [IN ]
  +                                     pertt,      &  ! [IN ]
  +                                     Xfact,      &  ! [IN ]
  +                                     lon,        &  ! [IN ]
  +                                     lat,        &  ! [IN ]
  +                                     p,          &  ! [IN ]
  +                                     z(k),       &  ! [IN ] (z)
  +                                     zcoords,    &  ! [IN ]
                                        wix(k),     &  ! [OUT] Zonal wind (m s^-1)(u)
                                        wiy(k),     &  ! [OUT] Meridional wind (m s^-1) (v)
                                        tmp(k),     &  ! [OUT] Temperature (K) (t)
                                        thetav(k),  &  ! [OUT] Virtual potential temperature (K)
  -                                     phis,    &  ! [OUT] Surface Geopotential (m^2 s^-2)
  -                                     ps,      &  ! [OUT] Surface Pressure (Pa)
  +                                     phis,       &  ! [OUT] Surface Geopotential (m^2 s^-2)
  +                                     ps,         &  ! [OUT] Surface Pressure (Pa)
                                        rho(k),     &  ! [OUT] density (kg m^-3)
                                        q(k)        )  ! [OUT] water vapor mixing ratio (kg/kg)
  -
  -       enddo
  -
  -!       ! iteration
  -!       do itr = 1, itrmax
  -
  -!          if ( itr == 1 ) then
  -!             eta(:,:) = 1.E-7_RP ! This initial value is recommended by Jablonowsky.
  -!          else
  -!             call eta_vert_coord_NW( kdim, itr, z_local, tmp, geo, eta_limit, eta, signal )
  -!          endif
  -
  -!          call steady_state( kdim, lat, eta, wix, wiy, tmp, geo )
  -
  -!          if( .NOT. signal ) exit
  +!          tv(k) = tmp(k) * (1.d0 + 0.61d0 * q(k))
  +          if ( z(k) < 1.0D-7 ) then
  +             prs(k) = ps
  +          else
  +             prs(k) = rho(k) * Rd * tmp(k)
  +          endif
  +       enddo
  +
  +!       ! pressure (upward: trapezoid)
  +!       prs(1) = ps
  +!       do k=2, kdim
  +!          dz = (z(k) - z(k-1))
  +!          prs(k) = prs(k-1) * ( 1.0_RP + dz*f_cf/(2.0_RP*Rd*tv(k-1)) ) &
  +!                            / ( 1.0_RP - dz*f_cf/(2.0_RP*Rd*tv(k)) )
   !       enddo
   
  -!       if ( itr > itrmax ) then
  -!          write(*          ,*) 'ETA ITERATION ERROR: NOT CONVERGED', n, l
  -!          write(ADM_LOG_FID,*) 'ETA ITERATION ERROR: NOT CONVERGED', n, l
  -!          call ADM_proc_stop
  -!       endif
  -
          call conv_vxvyvz ( kdim, lat, lon, wix, wiy, vx_local, vy_local, vz_local )
   
          do k = 1, kdim
  -          DIAG_var(n,k,l,1   ) = prs(k)
  -          DIAG_var(n,k,l,2   ) = tmp(k)
  -          DIAG_var(n,k,l,3   ) = vx_local(k)
  -          DIAG_var(n,k,l,4   ) = vy_local(k)
  -          DIAG_var(n,k,l,5   ) = vz_local(k)
  -          DIAG_var(n,k,l,I_QV) = q(k)
  +          DIAG_var(n,k,l,1     ) = real(prs(k),kind=RP)
  +          DIAG_var(n,k,l,2     ) = real(tmp(k),kind=RP)
  +          DIAG_var(n,k,l,3     ) = real(vx_local(k),kind=RP)
  +          DIAG_var(n,k,l,4     ) = real(vy_local(k),kind=RP)
  +          DIAG_var(n,k,l,5     ) = real(vz_local(k),kind=RP)
  +          DIAG_var(n,k,l,6+I_QV) = real(q(k),kind=RP)
          enddo
   
          if ( chemtracer ) then
  
positive.code:     integer,          intent(in)  :: lall
      character(len=*), intent(in)  :: test_case
  !    real(RP),         intent(in)  :: eps_geo2prs
  !    logical,          intent(in)  :: nicamcore
      logical,          intent(in)  :: chemtracer
      real(RP),         intent(out) :: DIAG_var(ijdim,kdim,lall,6+TRC_VMAX)
  
  
negative.code: 
            ! Todo : the mixing ratios are dry
            ! i.e. the ratio between the density of the species and the density of dry air.
            ! calc density at reference level
            DIAG_var(n,:,l,6+NCHEM_STR) = cl
            DIAG_var(n,:,l,6+NCHEM_END) = cl2
         endif
  
      enddo
      enddo
  
      return
    end subroutine jbw_moist_init
  