! careful adding any 1 stuff in here, adjust_tempqv has
! c3 and c4 to compute pressure with two reference pressures

SUBROUTINE init_domain_constants_em ( parent , nest )
   USE module_domain, ONLY : domain
   IMPLICIT NONE
   TYPE(domain)  :: parent , nest

   INTEGER iswater, islake, isice, isurban, isoilwater, map_proj, julyr, julday
   REAL    truelat1 , truelat2 , gmt , moad_cen_lat , stand_lon, pole_lat, pole_lon
   CHARACTER (LEN=256) :: char_junk

! single-value constants

   nest%p_top   = parent%p_top
   nest%save_topo_from_real   = parent%save_topo_from_real
   nest%cfn     = parent%cfn
   nest%cfn1    = parent%cfn1
   nest%rdx     = 1./nest%dx
   nest%rdy     = 1./nest%dy
!  nest%dts     = nest%dt/float(nest%time_step_sound)
   nest%dtseps  = parent%dtseps  ! used in height model only?
   nest%resm    = parent%resm    ! used in height model only?
   nest%zetatop = parent%zetatop ! used in height model only?
   nest%cf1     = parent%cf1
   nest%cf2     = parent%cf2
   nest%cf3     = parent%cf3
   nest%gmt     = parent%gmt
   nest%julyr   = parent%julyr
   nest%julday  = parent%julday
   nest%iswater = parent%iswater
   nest%isice   = parent%isice
   nest%isurban = parent%isurban
   nest%islake  = parent%islake
   nest%isoilwater = parent%isoilwater
   nest%mminlu  = trim(parent%mminlu)
   nest%tiso    = parent%tiso
   nest%tlp     = parent%tlp
   nest%p00     = parent%p00
   nest%t00     = parent%t00
   nest%tlp_strat= parent%tlp_strat
   nest%p_strat = parent%p_strat
!cyl: variables for trajectory /float
   nest%traj_k = parent%traj_k
   nest%traj_long = parent%traj_long
   nest%traj_lat = parent%traj_lat
   nest%this_is_an_ideal_run = parent%this_is_an_ideal_run
   nest%lake_depth_flag = parent%lake_depth_flag

   CALL nl_get_mminlu ( 1, char_junk )
   CALL nl_get_iswater( 1, iswater )
   CALL nl_get_islake ( 1, islake )
   CALL nl_get_isice  ( 1, isice )
   CALL nl_get_isurban( 1, isurban )
   CALL nl_get_isoilwater(1, isoilwater )
   CALL nl_get_truelat1 ( 1 , truelat1 )
   CALL nl_get_truelat2 ( 1 , truelat2 )
   CALL nl_get_moad_cen_lat ( 1 , moad_cen_lat )
   CALL nl_get_stand_lon ( 1 , stand_lon )
   CALL nl_get_pole_lat ( 1 , pole_lat )
   CALL nl_get_pole_lon ( 1 , pole_lon )
   CALL nl_get_map_proj ( 1 , map_proj )
   CALL nl_get_gmt ( 1 , gmt)
   CALL nl_get_julyr ( 1 , julyr)
   CALL nl_get_julday ( 1 , julday)
   IF ( nest%id .NE. 1 ) THEN
     CALL nl_set_gmt (nest%id, gmt)
     CALL nl_set_julyr (nest%id, julyr)
     CALL nl_set_julday (nest%id, julday)
     CALL nl_set_iswater ( nest%id, iswater )
     CALL nl_set_islake  ( nest%id, islake )
     CALL nl_set_isice   ( nest%id, isice )
     CALL nl_set_isurban ( nest%id, isurban )
     CALL nl_set_isoilwater ( nest%id, isoilwater )
     CALL nl_set_mminlu ( nest%id, char_junk )
     CALL nl_set_truelat1 ( nest%id , truelat1 )
     CALL nl_set_truelat2 ( nest%id , truelat2 )
     CALL nl_set_moad_cen_lat ( nest%id , moad_cen_lat )
     CALL nl_set_stand_lon ( nest%id , stand_lon )
     CALL nl_set_pole_lat ( nest%id , pole_lat )
     CALL nl_set_pole_lon ( nest%id , pole_lon )
     CALL nl_set_map_proj ( nest%id , map_proj )
   END IF
   nest%gmt     = gmt
   nest%julday  = julday
   nest%julyr   = julyr
   nest%iswater = iswater
   nest%islake  = islake
   nest%isice   = isice
   nest%isoilwater = isoilwater
   nest%mminlu  = trim(char_junk)
   nest%truelat1= truelat1
   nest%truelat2= truelat2
   nest%moad_cen_lat= moad_cen_lat
   nest%stand_lon= stand_lon
   nest%pole_lat= pole_lat
   nest%pole_lon= pole_lon
   nest%map_proj= map_proj

   nest%step_number  = parent%step_number

! 1D constants (Z)

!DAVE - maybe test if vert_nest instead
!  IF  ( parent%e_vert .EQ. nest%e_vert ) THEN
      nest%fnm(1:parent%e_vert)       = parent%fnm(1:parent%e_vert)
      nest%fnp(1:parent%e_vert)       = parent%fnp(1:parent%e_vert)
      nest%rdnw(1:parent%e_vert)      = parent%rdnw(1:parent%e_vert)
      nest%rdn(1:parent%e_vert)       = parent%rdn(1:parent%e_vert)
      nest%dnw(1:parent%e_vert)       = parent%dnw(1:parent%e_vert)
      nest%dn(1:parent%e_vert)        = parent%dn(1:parent%e_vert)
      nest%znu(1:parent%e_vert)       = parent%znu(1:parent%e_vert)
      nest%znw(1:parent%e_vert)       = parent%znw(1:parent%e_vert)
      nest%t_base(1:parent%e_vert)    = parent%t_base(1:parent%e_vert)
      nest%u_base(1:parent%e_vert)    = parent%u_base(1:parent%e_vert)
      nest%v_base(1:parent%e_vert)    = parent%v_base(1:parent%e_vert)
      nest%qv_base(1:parent%e_vert)   = parent%qv_base(1:parent%e_vert)
      nest%z_base(1:parent%e_vert)    = parent%z_base(1:parent%e_vert)
      nest%c1h(1:parent%e_vert)       = parent%c1h(1:parent%e_vert)
      nest%c2h(1:parent%e_vert)       = parent%c2h(1:parent%e_vert)
      nest%c3h(1:parent%e_vert)       = parent%c3h(1:parent%e_vert)
      nest%c4h(1:parent%e_vert)       = parent%c4h(1:parent%e_vert)
      nest%c1f(1:parent%e_vert)       = parent%c1f(1:parent%e_vert)
      nest%c2f(1:parent%e_vert)       = parent%c2f(1:parent%e_vert)
      nest%c3f(1:parent%e_vert)       = parent%c3f(1:parent%e_vert)
      nest%c4f(1:parent%e_vert)       = parent%c4f(1:parent%e_vert)
!  END IF

   nest%dzs       = parent%dzs
   nest%zs        = parent%zs

END SUBROUTINE init_domain_constants_em


!---------------------------------------------------------------------------------------------------

SUBROUTINE init_domain_vert_nesting ( parent, nest, use_baseparam_fr_nml )

!KAL this is a driver to initialize the vertical coordinates for the nest when vertical nesting is used
   USE module_domain
   IMPLICIT NONE
   TYPE(domain), POINTER                      :: parent, nest
   LOGICAL                                    :: use_baseparam_fr_nml
   !local
   REAL, DIMENSION(parent%e_vert)             :: znw_c

   INTERFACE

      SUBROUTINE vert_cor_vertical_nesting_integer(nest,znw_c,k_dim_c)
         USE module_domain
         TYPE(domain), POINTER :: nest
         integer , intent(in) :: k_dim_c
         real , dimension(k_dim_c), INTENT(IN) :: znw_c
      END SUBROUTINE vert_cor_vertical_nesting_integer

      SUBROUTINE vert_cor_vertical_nesting_arbitrary(nest,znw_c,kde_c,use_baseparam_fr_nml)
          USE module_domain
          TYPE(domain), POINTER :: nest
          INTEGER, INTENT(IN   ) :: kde_c
          REAL, DIMENSION(kde_c), INTENT(IN   ) :: znw_c
          LOGICAL, INTENT(IN   ) :: use_baseparam_fr_nml
      END SUBROUTINE vert_cor_vertical_nesting_arbitrary
      
   END INTERFACE

   
   ! save the coarse grid values here
   
   znw_c = nest%znw(1:parent%e_vert)
   
   ! calculate the nest (fine) grid values here
   ! one of these calls goes to integer refinement in the vertical direction, and one goes to arbitrary refinement.  Eventually the call to integer refinement will be obsolete.
   if (nest%vert_refine_method .EQ. 1) then  !if you are in this subroutine there is vertical nesting- (i.e. nest%e_vert /= parent%e_vert to enter this subroutine)
      CALL vert_cor_vertical_nesting_integer(nest,znw_c,parent%e_vert)
   elseif (nest%vert_refine_method .EQ. 2) then
      CALL vert_cor_vertical_nesting_arbitrary(nest,znw_c,parent%e_vert,use_baseparam_fr_nml)
   endif

END SUBROUTINE init_domain_vert_nesting

!-----------------------------------------------------------------------------------------


!this is a direct copy of a subroutine that is in ndown, but I couldn't link to the subroutine in ndown because it is compiled after this file
!so a dependecy on ndown will not work. Additionally, ndown is not compiled for ideal cases. The variable is named parent in ndown, but it is actually operating on the nest. It has been renamed to nest here.
      SUBROUTINE vert_cor_vertical_nesting_integer(nest,znw_c,k_dim_c)
         USE module_domain
         USE module_model_constants
   IMPLICIT NONE
         TYPE(domain), POINTER :: nest
         integer , intent(in) :: k_dim_c
         real , dimension(k_dim_c), INTENT(IN) :: znw_c

       integer :: kde_c , kde_n ,n_refine,ii,kkk,k
       real :: dznw_m,cof1,cof2

       INTEGER :: ids, ide, jds, jde, kds, kde, &
                  ims, ime, jms, jme, kms, kme, &
                  its, ite, jts, jte, kts, kte

!KAL this subroutine recalculates the vertical coordinates for the nest when vertical nesting is used.  This routine is copied from ndown and allows integer refinement only.

!KAL znw is eta values on full w levels
!KAL everything else is set from znw
!KAL dnw is delta eta on w levels
!KAL rdn is inverse delta eta on w levels
!KAL fnp

        kde_c = k_dim_c
        kde_n = nest%e_vert
!        n_refine = nest%vert_refine_fact
        n_refine = (kde_n-1)/(kde_c-1)

         kkk = 0
         do k = 1 , kde_c-1
         dznw_m = znw_c(k+1) - znw_c(k)
         do ii = 1,n_refine
         kkk = kkk + 1
         nest%znw(kkk) = znw_c(k) + float(ii-1)/float(n_refine)*dznw_m
         enddo
         enddo
         nest%znw(kde_n) = znw_c(kde_c)
         nest%znw(1) = znw_c(1)

      DO k=1, kde_n-1
         nest%dnw(k) = nest%znw(k+1) - nest%znw(k)
         nest%rdnw(k) = 1./nest%dnw(k)
         nest%znu(k) = 0.5*(nest%znw(k+1)+nest%znw(k))
      END DO

      DO k=2, kde_n-1
         nest%dn(k) = 0.5*(nest%dnw(k)+nest%dnw(k-1))
         nest%rdn(k) = 1./nest%dn(k)
         nest%fnp(k) = .5* nest%dnw(k  )/nest%dn(k)
         nest%fnm(k) = .5* nest%dnw(k-1)/nest%dn(k)
      END DO

  cof1 = (2.*nest%dn(2)+nest%dn(3))/(nest%dn(2)+nest%dn(3))*nest%dnw(1)/nest%dn(2)
  cof2 =     nest%dn(2)        /(nest%dn(2)+nest%dn(3))*nest%dnw(1)/nest%dn(3)

      nest%cf1  = nest%fnp(2) + cof1
      nest%cf2  = nest%fnm(2) - cof1 - cof2
      nest%cf3  = cof2

      nest%cfn  = (.5*nest%dnw(kde_n-1)+nest%dn(kde_n-1))/nest%dn(kde_n-1)
      nest%cfn1 = -.5*nest%dnw(kde_n-1)/nest%dn(kde_n-1)
      
      ! the variables dzs and zs are kept from the parent domain.  These are the depths and thickness of the soil layers, which are not included in vertical nesting.

    CALL get_ijk_from_grid( nest, &
                            ids, ide, jds, jde, kds, kde, &
                            ims, ime, jms, jme, kms, kme, &
                            its, ite, jts, jte, kts, kte )

    CALL compute_vcoord_1d_coeffs ( nest%ht, nest%etac, nest%znw, &
                                    model_config_rec%hybrid_opt, &
                                    r_d, g, p1000mb, &
                                    nest%p_top, nest%p00, nest%t00, nest%tlp, &
                                    ids, ide, jds, jde, kds, kde, &
                                    ims, ime, jms, jme, kms, kme, &
                                    its, ite, jts, jte, kts, kte, &
                                    nest%znu, &
                                    nest%c1f, nest%c2f, nest%c3f, nest%c4f, &
                                    nest%c1h, nest%c2h, nest%c3h, nest%c4h )

      END SUBROUTINE vert_cor_vertical_nesting_integer

!-----------------------------------------------------------------------------------------

SUBROUTINE vert_cor_vertical_nesting_arbitrary(nest,znw_c,kde_c,use_baseparam_fr_nml)
    USE module_domain
    USE module_model_constants
    IMPLICIT NONE
    TYPE(domain), POINTER :: nest
    INTEGER, INTENT(IN   ) :: kde_c
    REAL, DIMENSION(kde_c), INTENT(IN   ) :: znw_c
    LOGICAL, INTENT(IN   ) :: use_baseparam_fr_nml
    INTEGER :: k, kde_n, ks, id
    REAL :: cof1, cof2
    REAL :: max_dz = 1000
    REAL :: p00, t00, a, tiso, a_strat, p_strat
    INTEGER :: ids, ide, jds, jde, kds, kde, &
               ims, ime, jms, jme, kms, kme, &
               its, ite, jts, jte, kts, kte
    REAL, DIMENSION(max_eta) :: eta_levels

    IF ( use_baseparam_fr_nml ) then
      CALL nl_get_base_pres        ( 1 , p00        )
      CALL nl_get_base_temp        ( 1 , t00        )
      CALL nl_get_base_lapse       ( 1 , a          )
      CALL nl_get_iso_temp         ( 1 , tiso       )
      CALL nl_get_base_lapse_strat ( 1 , a_strat    )
      CALL nl_get_base_pres_strat  ( 1 , p_strat    )
      IF ((t00 .LT. 100.0) .OR. (p00 .LT. 10000.0)) THEN
        WRITE(wrf_err_message,*) '--- ERROR: bad base state for T00 or P00 in namelist.input file'
        CALL wrf_error_fatal(TRIM(wrf_err_message))
      END IF
    ELSE
      p00     = nest%p00
      t00     = nest%t00
      a       = nest%tlp
      tiso    = nest%tiso
      a_strat = nest%tlp_strat
      p_strat = nest%p_strat
      IF ((t00 .LT. 100.0) .OR. (p00 .LT. 10000.0)) THEN
        WRITE(wrf_err_message,*) '--- ERROR: did not find base state parameters in nest. Add use_baseparam_fr_nml = .true. in &dynamics and rerun'
        CALL wrf_error_fatal(TRIM(wrf_err_message))
      ENDIF
    ENDIF

    kde_n = nest%e_vert

    !DJW Added code for specifying multiple domains' eta_levels
    IF (nest%id .NE. 1) THEN
      id = 1
      ks = 1
      DO WHILE (nest%id .GT. id)
        id = id+1
        ks = ks+model_config_rec%e_vert(id-1)
      ENDDO
    ENDIF
    IF ((nest%this_is_an_ideal_run) .AND. (model_config_rec%eta_levels(1) .EQ. -1.0)) THEN
      !DJW If we're running an ideal case and do not have levels set in the
      !namelist then we set znw using constant spacing in eta
      DO k=1,kde_n
        CALL wrf_debug(0, "nest_init_utils: eta_levels are not specified in the namelist, setting levels with constant spacing in eta.")
        nest%znw(k) = 1.0-(k-1)/FLOAT((kde_n-1))
      ENDDO
    ELSEIF (.NOT.(nest%this_is_an_ideal_run) .AND. (model_config_rec%eta_levels(1) .EQ. -1.0)) THEN
      write(*,'(A,I2,A)') "--- WARNING: eta_levels are not specified in the namelist for grid_id=",nest%grid_id,", using WRF's default levels."
      CALL get_ijk_from_grid( nest, &
                              ids, ide, jds, jde, kds, kde, &
                              ims, ime, jms, jme, kms, kme, &
                              its, ite, jts, jte, kts, kte )
      write(*,'(A,F10.3)') "--- USING: nest%p_top   = ",nest%p_top
      write(*,'(A,F10.3)') "--- USING: g            = ",g
      write(*,'(A,F10.3)') "--- USING: cvpm         = ",cvpm
      write(*,'(A,F10.3)') "--- USING: r_d          = ",r_d
      write(*,'(A,F10.3)') "--- USING: cp           = ",cp
      write(*,'(A,F10.3)') "--- USING: p1000mb      = ",p1000mb
      write(*,'(A,F10.3)') "--- USING: t0           = ",t0
      write(*,'(A,F10.3)') "--- USING: p00          = ",p00
      write(*,'(A,F10.3)') "--- USING: t00          = ",t00
      write(*,'(A,F10.3)') "--- USING: a            = ",a
      write(*,'(A,F10.3)') "--- USING: tiso         = ",tiso
      write(*,'(A,F10.3)') "--- USING: a_strat      = ",a_strat
      write(*,'(A,F10.3)') "--- USING: p_strat      = ",p_strat
      CALL compute_eta ( nest%znw, &
                         eta_levels, max_eta, max_dz, &
                         nest%p_top, g, p00, cvpm, a, r_d, cp, &
                         t00, p1000mb, t0, tiso, p_strat, a_strat, &
                         ids, ide, jds, jde, kds, kde, &
                         ims, ime, jms, jme, kms, kme, &
                         its, ite, jts, jte, kts, kte )
    ELSE
      !DJW If we're in here then we are suppose to read in eta_levels
      !from the namelist. We do so and then check to make sure they make sense
    DO k=1,kde_n
      nest%znw(k) = model_config_rec%eta_levels(ks+k-1)
        write(*,'(A,I3,A,F6.3)') "nest%znw(",k,") = ",nest%znw(k)
    ENDDO
    !Check the value of the first and last eta level for our domain,
    !then check that the vector of eta levels is only decreasing
    IF (nest%znw(1) .NE. 1.0) THEN
        write(wrf_err_message,'(A,I2,A)') "--- ERROR: first eta_level for grid_id=",nest%grid_id," is not 1.0. Check namelist."
        CALL wrf_error_fatal( wrf_err_message )
    ENDIF
      write(*,'(A,F10.3)') "--- USING: g            = ",g
      write(*,'(A,F10.3)') "--- USING: cvpm         = ",cvpm
      write(*,'(A,F10.3)') "--- USING: r_d          = ",r_d
      write(*,'(A,F10.3)') "--- USING: cp           = ",cp
      write(*,'(A,F10.3)') "--- USING: p1000mb      = ",p1000mb
      write(*,'(A,F10.3)') "--- USING: t0           = ",t0
    IF (nest%znw(kde_n) .NE. 0.0) THEN
        write(wrf_err_message,'(A,I2,A)') "--- ERROR: last eta_level for grid_id=",nest%grid_id," is not 0.0. Check namelist."
        CALL wrf_error_fatal( wrf_err_message )
    ENDIF
    DO k=2,kde_n
      IF (nest%znw(k) .GT. nest%znw(k-1)) THEN
          write(wrf_err_message,'(A,I2,A)') "--- ERROR: eta_level for grid_id=",nest%grid_id," are not monotonically decreasing. Check namelist."
          CALL wrf_error_fatal( wrf_err_message )
      ENDIF
    ENDDO
    ENDIF
    !DJW End of added code for specifying eta_levels

    DO k=1,kde_n-1
      nest%dnw(k) = nest%znw(k+1)-nest%znw(k)
      nest%rdnw(k) = 1./nest%dnw(k)
      nest%znu(k) = 0.5*(nest%znw(k+1)+nest%znw(k))
    ENDDO
    nest%znu(kde_n) = 0.0

    DO k=2,kde_n-1
      nest%dn(k) = 0.5*(nest%dnw(k)+nest%dnw(k-1))
      nest%rdn(k) = 1./nest%dn(k)
      nest%fnp(k) = .5* nest%dnw(k  )/nest%dn(k)
      nest%fnm(k) = .5* nest%dnw(k-1)/nest%dn(k)
    ENDDO

    cof1 = (2.*nest%dn(2)+nest%dn(3))/(nest%dn(2)+nest%dn(3))*nest%dnw(1)/nest%dn(2)
    cof2 =     nest%dn(2)            /(nest%dn(2)+nest%dn(3))*nest%dnw(1)/nest%dn(3)

    nest%cf1  = nest%fnp(2) + cof1
    nest%cf2  = nest%fnm(2) - cof1 - cof2
    nest%cf3  = cof2
    nest%cfn  = (.5*nest%dnw(kde_n-1)+nest%dn(kde_n-1))/nest%dn(kde_n-1)
    nest%cfn1 = -.5*nest%dnw(kde_n-1)/nest%dn(kde_n-1)

    CALL get_ijk_from_grid( nest, &
                            ids, ide, jds, jde, kds, kde, &
                            ims, ime, jms, jme, kms, kme, &
                            its, ite, jts, jte, kts, kte )

    CALL compute_vcoord_1d_coeffs ( nest%ht, nest%etac, nest%znw, &
                                    model_config_rec%hybrid_opt, &
                                    r_d, g, p1000mb, &
                                    nest%p_top, nest%p00, nest%t00, nest%tlp, &
                                    ids, ide, jds, jde, kds, kde, &
                                    ims, ime, jms, jme, kms, kme, &
                                    its, ite, jts, jte, kts, kte, &
                                    nest%znu, &
                                    nest%c1f, nest%c2f, nest%c3f, nest%c4f, &
                                    nest%c1h, nest%c2h, nest%c3h, nest%c4h )
    
END SUBROUTINE vert_cor_vertical_nesting_arbitrary

SUBROUTINE compute_eta ( znw , &
                         eta_levels , max_eta , max_dz ,        &
                         p_top_def , g_def , p00_def ,          &
                         cvpm_def , a_def , r_d_def , cp_def ,  &
                         t00_def , p1000mb_def , t0_def ,       &
                         tiso_def , p_strat_def , a_strat_def , &
                         ids , ide , jds , jde , kds , kde ,    &
                         ims , ime , jms , jme , kms , kme ,    &
                         its , ite , jts , jte , kts , kte )

      !  Compute eta levels, either using given values from the namelist (hardly
      !  a computation, yep, I know), or assuming a constant dz above the PBL,
      !  knowing p_top and the number of eta levels.

      IMPLICIT NONE

      INTEGER , INTENT(IN)        :: ids , ide , jds , jde , kds , kde , &
                                     ims , ime , jms , jme , kms , kme , &
                                     its , ite , jts , jte , kts , kte
      REAL , INTENT(IN)           :: max_dz
      REAL , INTENT(IN)           :: p_top_def , g_def , p00_def , cvpm_def , &
                                     a_def , r_d_def , cp_def , t00_def ,     &
                                     p1000mb_def , t0_def , tiso_def ,        &
                                     p_strat_def , a_strat_def
      INTEGER , INTENT(IN)        :: max_eta
      REAL , DIMENSION (max_eta)  :: eta_levels

      REAL , DIMENSION (kts:kte) , INTENT(OUT) :: znw

      !  Local vars

      INTEGER :: k , kk
      REAL(KIND=8) :: mub , t_init , p_surf , pb, ztop, ztop_pbl , dz , temp
      REAL(KIND=8) , DIMENSION(kts:kte) :: dnw
      REAL(KIND=8) :: p_top , g , p00 , cvpm , &
                      a , r_d , cp , t00 ,     &
                      p1000mb , t0 , tiso ,        &
                      p_strat , a_strat

      INTEGER , PARAMETER :: prac_levels = 59
      INTEGER :: loop , loop1
      REAL(KIND=8) , DIMENSION(prac_levels) :: znw_prac , znu_prac , dnw_prac
      REAL(KIND=8) , DIMENSION(MAX(prac_levels,kde)) :: alb , phb
      REAL(KIND=8) :: alb_max, t_init_max, pb_max, phb_max

      CHARACTER(LEN=256) :: message

      !  Compute top of the atmosphere with some silly levels.  We just want to
      !  integrate to get a reasonable value for ztop.  We use the planned
      !  PBL-esque
      !  levels, and then just coarse resolution above that.  We know p_top,
      !  and we
      !  have the base state vars.

      p_top   = p_top_def
      g       = g_def
      p00     = p00_def
      cvpm    = cvpm_def
      a       = a_def
      r_d     = r_d_def
      cp      = cp_def
      t00     = t00_def
      p1000mb = p1000mb_def
      t0      = t0_def
      tiso    = tiso_def
      p_strat = p_strat_def
      a_strat = a_strat_def

      p_surf = p00

      znw_prac = (/ 1.0000_8 , 0.9930_8 , 0.9830_8 , 0.9700_8 , 0.9540_8 , 0.9340_8 , 0.9090_8 , 0.8800_8 , &
                    0.8500_8 , 0.8000_8 , 0.7500_8 , 0.7000_8 , 0.6500_8 , 0.6000_8 , 0.5500_8 , 0.5000_8 , &
                    0.4500_8 , 0.4000_8 , 0.3500_8 , 0.3000_8 , 0.2500_8 , 0.2000_8 , 0.1500_8 , 0.1000_8 , &
                    0.0800_8 , 0.0600_8 , 0.0400_8 , 0.0200_8 , &
                    0.0150_8 , 0.0100_8 , 0.0090_8 , 0.0080_8 , 0.0070_8 , 0.0060_8 , 0.0050_8 , 0.0040_8 , &
                    0.0035_8 , 0.0030_8 , &
                    0.0028_8 , 0.0026_8 , 0.0024_8 , 0.0022_8 , 0.0020_8 , &
                    0.0018_8 , 0.0016_8 , 0.0014_8 , 0.0012_8 , 0.0010_8 , &
                    0.0009_8 , 0.0008_8 , 0.0007_8 , 0.0006_8 , 0.0005_8 , 0.0004_8 , 0.0003_8 , &
                    0.0002_8 , 0.0001_8 , 0.00005_8, 0.0000_8 /)

      DO k = 1 , prac_levels - 1
        znu_prac(k) = ( znw_prac(k) + znw_prac(k+1) ) * 0.5_8
        dnw_prac(k) = znw_prac(k+1) - znw_prac(k)
      END DO

      DO k = 1, prac_levels-1
        pb = znu_prac(k)*(p_surf - p_top) + p_top
        temp = MAX ( tiso, t00 + A*LOG(pb/p00) )
        IF ( pb .LT. p_strat ) THEN
          temp = tiso + A_strat*LOG(pb/p_strat)
        END IF
        t_init = temp*(p00/pb)**(r_d/cp) - t0
        alb(k) = (r_d/p1000mb)*(t_init+t0)*(pb/p1000mb)**cvpm
      END DO

      !  Base state mu is defined as base state surface pressure minus p_top

      mub = p_surf - p_top

      !  Integrate base geopotential, starting at terrain elevation.

      phb(1) = 0._8
      DO k  = 2,prac_levels
        phb(k) = phb(k-1) - dnw_prac(k-1)*mub*alb(k-1)
      END DO

      !  So, now we know the model top in meters.  Get the average depth above the PBL
      !  of each of the remaining levels.  We are going for a constant delta z thickness.

      ztop     = phb(prac_levels) / g
      ztop_pbl = phb(8          ) / g
      dz = ( ztop - ztop_pbl ) / REAL ( kde - 8 )

      IF ( dz .GE. max_dz ) THEN
        WRITE (message,FMT='("With a requested ",F7.1," Pa model top, the model lid will be about ",F7.1," m.")') p_top, ztop
        CALL wrf_message ( message )
        WRITE (message,FMT='("With ",I3," levels above the PBL, the level thickness will be about ",F6.1," m.")') kde-8, dz
        CALL wrf_message ( message )
        WRITE (message,FMT='("Thicknesses greater than ",F7.1," m are not recommended.")') max_dz
        CALL wrf_message ( message )
        CALL wrf_error_fatal ( 'Add more levels to namelist.input for e_vert' )
      END IF

      !  Standard levels near the surface so no one gets in trouble.

      DO k = 1 , 8
        eta_levels(k) = znw_prac(k)
      END DO

      !  Using d phb(k)/ d eta(k) = -mub * alb(k), eqn 2.9
      !  Skamarock et al, NCAR TN 468.  Use full levels, so
      !  use twice the thickness.

      DO k = 8, kte-1-2

        find_prac : DO kk = 1 , prac_levels
          IF (znw_prac(kk) .LT. eta_levels(k) ) THEN
            EXIT find_prac
          END IF
        end do find_prac

        pb = 0.5*(eta_levels(k)+znw_prac(kk)) * (p_surf - p_top) + p_top

        temp = MAX ( tiso, t00 + A*LOG(pb/p00) )
        IF ( pb .LT. p_strat ) THEN
          temp = tiso + A_strat * LOG ( pb/p_strat )
        END IF
!       temp =             t00 + A*LOG(pb/p00)
        t_init = temp*(p00/pb)**(r_d/cp) - t0
        alb(k) = (r_d/p1000mb)*(t_init+t0)*(pb/p1000mb)**cvpm
        eta_levels(k+1) = eta_levels(k) - dz*g / ( mub*alb(k) )
        pb = 0.5*(eta_levels(k)+eta_levels(k+1)) * (p_surf - p_top) + p_top

        temp = MAX ( tiso, t00 + A*LOG(pb/p00) )
        IF ( pb .LT. p_strat ) THEN
          temp = tiso + A_strat * LOG ( pb/p_strat )
        END IF
!       temp =             t00 + A*LOG(pb/p00)
        t_init = temp*(p00/pb)**(r_d/cp) - t0
        alb(k) = (r_d/p1000mb)*(t_init+t0)*(pb/p1000mb)**cvpm
        eta_levels(k+1) = eta_levels(k) - dz*g / ( mub*alb(k) )
        pb = 0.5*(eta_levels(k)+eta_levels(k+1)) * (p_surf - p_top) + p_top

        phb(k+1) = phb(k) - (eta_levels(k+1)-eta_levels(k)) * mub*alb(k)
      END DO

      alb_max = alb(kte-1-2)
      t_init_max = t_init
      pb_max = pb
      phb_max = phb(kte-1)

      DO k = 1 , kte-1-2
        znw(k) = eta_levels(k)
      END DO
      znw(kte-2) = 0.000

      !  There is some iteration.  We want the top level, ztop, to be
      !  consistent with the delta z, and we want the half level values
      !  to be consistent with the eta levels.  The inner loop to 10 gets
      !  the eta levels very accurately, but has a residual at the top, due
      !  to dz changing.  We reset dz five times, and then things seem OK.

      DO loop1 = 1 , 5
        DO loop = 1 , 10
          DO k = 8, kte-1-2-1
            pb = (znw(k)+znw(k+1))*0.5 * (p_surf - p_top) + p_top
            temp = MAX ( tiso, t00 + A*LOG(pb/p00) )
            IF ( pb .LT. p_strat ) THEN
              temp = tiso + A_strat * LOG ( pb/p_strat )
            END IF
            t_init = temp*(p00/pb)**(r_d/cp) - t0
            alb(k) = (r_d/p1000mb)*(t_init+t0)*(pb/p1000mb)**cvpm
            znw(k+1) = znw(k) - dz*g / ( mub*alb(k) )
          END DO
          pb = pb_max
          t_init = t_init_max
          alb(kte-1-2) = alb_max
          znw(kte-2) = znw(kte-1-2) - dz*g / ( mub*alb(kte-1-2) )
          IF ( ( loop1 .EQ. 5 ) .AND. ( loop .EQ. 10 ) ) THEN
            print *,'Converged znw(kte) should be about 0.0 = ',znw(kte-2)
          END IF
          znw(kte-2) = 0.000
        END DO

        !  Here is where we check the eta levels values we just computed.

        DO k = 1, kde-1-2
          pb = (znw(k)+znw(k+1))*0.5 * (p_surf - p_top) + p_top
          temp = MAX ( tiso, t00 + A*LOG(pb/p00) )
          IF ( pb .LT. p_strat ) THEN
            temp = tiso + A_strat * LOG ( pb/p_strat )
          END IF
          t_init = temp*(p00/pb)**(r_d/cp) - t0
          alb(k) = (r_d/p1000mb)*(t_init+t0)*(pb/p1000mb)**cvpm
        END DO

        phb(1) = 0.
        DO k  = 2,kde-2
          phb(k) = phb(k-1) - (znw(k)-znw(k-1)) * mub*alb(k-1)
        END DO

        !  Reset the model top and the dz, and iterate.

        ztop = phb(kde-2)/g
        ztop_pbl = phb(8)/g
        dz = ( ztop - ztop_pbl ) / REAL ( (kde-2) - 8 )
      END DO

      IF ( dz .GT. max_dz ) THEN
        print *,'z (m)            = ',phb(1)/g
        do k = 2 ,kte-2
          print *,'z (m) and dz (m) = ',phb(k)/g,(phb(k)-phb(k-1))/g
        end do
        print *,'dz (m) above fixed eta levels = ',dz
        print *,'namelist max_dz (m) = ',max_dz
        print *,'namelist p_top (Pa) = ',p_top
        CALL wrf_debug ( 0, 'You need one of three things:' )
        CALL wrf_debug ( 0, '1) More eta levels to reduce the dz: e_vert' )
        CALL wrf_debug ( 0, '2) A lower p_top so your total height is reduced: p_top_requested')
        CALL wrf_debug ( 0, '3) Increase the maximum allowable eta thickness: max_dz')
        CALL wrf_debug ( 0, 'All are namelist options')
        CALL wrf_error_fatal ( 'dz above fixed eta levels is too large')
      END IF

      !  Add those 2 levels back into the middle, just above the 8 levels
      !  that semi define a boundary layer.  After we open up the levels,
      !  then we just linearly interpolate in znw.  So now levels 1-8 are
      !  specified as the fixed boundary layer levels given in this routine.
      !  The top levels, 12 through kte are those computed.  The middle
      !  levels 9, 10, and 11 are equi-spaced in znw, and are each 1/2 the
      !  the znw thickness of levels 11 through 12.

      DO k = kte-2 , 9 , -1
        znw(k+2) = znw(k)
      END DO

      znw( 9) = 0.75 * znw( 8) + 0.25 * znw(12)
      znw(10) = 0.50 * znw( 8) + 0.50 * znw(12)
      znw(11) = 0.25 * znw( 8) + 0.75 * znw(12)

      DO k = 8, kte-1
        pb = (znw(k)+znw(k+1))*0.5 * (p_surf - p_top) + p_top
        temp = MAX ( tiso, t00 + A*LOG(pb/p00) )
        IF ( pb .LT. p_strat ) THEN
          temp = tiso + A_strat * LOG ( pb/p_strat )
        END IF
        t_init = temp*(p00/pb)**(r_d/cp) - t0
        alb(k) = (r_d/p1000mb)*(t_init+t0)*(pb/p1000mb)**cvpm
        phb(k) = phb(k-1) - (znw(k)-znw(k-1)) * mub*alb(k-1)
      END DO
      phb(kte) = phb(kte-1) - (znw(kte)-znw(kte-1)) * mub*alb(kte-1)

      k=1
      WRITE (*,FMT='("Full level index = ",I4,"     Height = ",F7.1," m")') k,phb(1)/g
      do k = 2 ,kte
        WRITE (*,FMT='("Full level index = ",I4,"     Height = ",F7.1," m      Thickness = ",F6.1," m")') k,phb(k)/g,(phb(k)-phb(k-1))/g
      end do

END SUBROUTINE compute_eta

!-----------------------------------------------------------------------------------------

SUBROUTINE blend_terrain ( ter_interpolated , ter_input , &
                           ids , ide , jds , jde , kds , kde , &
                           ims , ime , jms , jme , kms , kme , &
                           ips , ipe , jps , jpe , kps , kpe )

   USE module_configure
   IMPLICIT NONE

   INTEGER , INTENT(IN)                       :: ids , ide , jds , jde , kds , kde , &
                                                 ims , ime , jms , jme , kms , kme , &
                                                 ips , ipe , jps , jpe , kps , kpe
   REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN)    :: ter_interpolated
   REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(INOUT) :: ter_input

   REAL , DIMENSION(ims:ime,kms:kme,jms:jme) :: ter_temp
   INTEGER :: i , j , k , spec_bdy_width
   REAL    :: r_blend_zones
   INTEGER blend_cell, blend_width

   !  The fine grid elevation comes from the horizontally interpolated
   !  parent elevation for the first spec_bdy_width row/columns, so we need
   !  to get that value.  We blend the coarse and fine in the next blend_width
   !  rows and columns.  After that, in the interior, it is 100% fine grid.

   CALL nl_get_spec_bdy_width ( 1, spec_bdy_width)
   CALL nl_get_blend_width ( 1, blend_width)

   !  Initialize temp values to the nest ter elevation.  This fills in the values
   !  that will not be modified below.

   DO j = jps , MIN(jpe, jde-1)
      DO k = kps , kpe
         DO i = ips , MIN(ipe, ide-1)
            ter_temp(i,k,j) = ter_input(i,k,j)
         END DO
      END DO
   END DO

   !  To avoid some tricky indexing, we fill in the values inside out.  This allows
   !  us to overwrite incorrect assignments.  There are replicated assignments, and
   !  there is much unnecessary "IF test inside of a loop" stuff.  For a large
   !  domain, this is only a patch; for a small domain, this is not a biggy.

   r_blend_zones = 1./(blend_width+1)
   DO j = jps , MIN(jpe, jde-1)
      DO k = kps , kpe
         DO i = ips , MIN(ipe, ide-1)
            DO blend_cell = blend_width,1,-1
               IF   ( ( i .EQ.       spec_bdy_width + blend_cell ) .OR.  ( j .EQ.       spec_bdy_width + blend_cell ) .OR. &
                      ( i .EQ. ide - spec_bdy_width - blend_cell ) .OR.  ( j .EQ. jde - spec_bdy_width - blend_cell ) ) THEN
                  ter_temp(i,k,j) = ( (blend_cell)*ter_input(i,k,j) + (blend_width+1-blend_cell)*ter_interpolated(i,k,j) ) &
                                    * r_blend_zones
               END IF
            ENDDO
            IF      ( ( i .LE.       spec_bdy_width     ) .OR.  ( j .LE.       spec_bdy_width     ) .OR. &
                      ( i .GE. ide - spec_bdy_width     ) .OR.  ( j .GE. jde - spec_bdy_width     ) ) THEN
               ter_temp(i,k,j) =      ter_interpolated(i,k,j)
            END IF
         END DO
      END DO
   END DO

   !  Set nest elevation with temp values.  All values not overwritten in the above
   !  loops have been previously set in the initial assignment.

   DO j = jps , MIN(jpe, jde-1)
      DO k = kps , kpe
         DO i = ips , MIN(ipe, ide-1)
            ter_input(i,k,j) = ter_temp(i,k,j)
         END DO
      END DO
   END DO

END SUBROUTINE blend_terrain

SUBROUTINE copy_3d_field ( ter_interpolated , ter_input , &
                           ids , ide , jds , jde , kds , kde , &
                           ims , ime , jms , jme , kms , kme , &
                           ips , ipe , jps , jpe , kps , kpe )

   IMPLICIT NONE

   INTEGER , INTENT(IN)                       :: ids , ide , jds , jde , kds , kde , &
                                                 ims , ime , jms , jme , kms , kme , &
                                                 ips , ipe , jps , jpe , kps , kpe
   REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(OUT) :: ter_interpolated
   REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN)  :: ter_input

   INTEGER :: i , j , k

   DO j = jps , MIN(jpe, jde-1)
      DO k = kps , kpe
         DO i = ips , MIN(ipe, ide-1)
            ter_interpolated(i,k,j) = ter_input(i,k,j)
         END DO
      END DO
   END DO

END SUBROUTINE copy_3d_field

SUBROUTINE adjust_tempqv ( mub, save_mub, c3, c4, znw, p_top, &
                           th, pp, qv,  &
                           use_theta_m, &
                           ids , ide , jds , jde , kds , kde , &
                           ims , ime , jms , jme , kms , kme , &
                           ips , ipe , jps , jpe , kps , kpe )

   !USE module_configure
   !USE module_domain
   USE module_model_constants

   !USE module_bc
   !USE module_io_domain
   !USE module_state_description
   !USE module_timing
   !USE module_soil_pre
   IMPLICIT NONE

   INTEGER , INTENT(IN)                       :: ids , ide , jds , jde , kds , kde , &
                                                 ims , ime , jms , jme , kms , kme , &
                                                 ips , ipe , jps , jpe , kps , kpe
   REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN)    :: mub, save_mub
   REAL , DIMENSION(kms:kme) , INTENT(IN)    :: znw
   REAL , DIMENSION(kms:kme) , INTENT(IN)    :: c3, c4
   REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(INOUT) :: th, pp, qv
   INTEGER , INTENT(IN)                       :: use_theta_m

   REAL , DIMENSION(ims:ime,kms:kme,jms:jme) :: p_old, p_new, rh
   REAL :: es,dth,tc,e,dth1,thloc
   INTEGER :: i , j , k

   real p_top


! p_old = full pressure before terrain blending; also compute initial RH
! which is going to be conserved during terrain blending
   DO j = jps , MIN(jpe, jde-1)
      DO k = kps , kpe-1
         DO i = ips , MIN(ipe, ide-1)
            p_old(i,k,j) = c4(k) +         c3(k)*save_mub(i,j) + p_top + pp(i,k,j)
            IF ( use_theta_m .EQ. 1 ) THEN
               tc = (th(i,k,j)+300.)*(p_old(i,k,j)/1.e5)**(2./7.)/(1.+R_v/R_d*qv(i,k,j)) - 273.15
            ELSE
               tc = (th(i,k,j)+300.)*(p_old(i,k,j)/1.e5)**(2./7.)                        - 273.15
            END IF
            es = 610.78*exp(17.0809*tc/(234.175+tc))
            e = qv(i,k,j)*p_old(i,k,j)/(0.622+qv(i,k,j))
            rh(i,k,j) = e/es
         END DO
      END DO
   END DO
! p_new = full pressure after terrain blending; also compute temperature correction and convert RH back to QV
   DO j = jps , MIN(jpe, jde-1)
      DO k = kps , kpe-1
         DO i = ips , MIN(ipe, ide-1)
            p_new(i,k,j) = c4(k) +         c3(k)*mub(i,j) + p_top + pp(i,k,j)
! 2*(g/cp-6.5e-3)*R_dry/g = -191.86e-3
            IF ( use_theta_m .EQ. 1 ) THEN
               thloc = (th(i,k,j)+300.)/(1.+R_v/R_d*qv(i,k,j))
            ELSE
               thloc = (th(i,k,j)+300.)
            END IF
            dth1 = -191.86e-3*(thloc         )/(p_new(i,k,j)+p_old(i,k,j))*(p_new(i,k,j)-p_old(i,k,j))
            dth  = -191.86e-3*(thloc+0.5*dth1)/(p_new(i,k,j)+p_old(i,k,j))*(p_new(i,k,j)-p_old(i,k,j))
            IF ( use_theta_m .EQ. 1 ) THEN
               th(i,k,j) = (thloc+dth)*(1.+R_v/R_d*qv(i,k,j))-300.
            ELSE
               th(i,k,j) = (thloc+dth)                       -300.
            END IF
            tc = (thloc+dth)*(p_new(i,k,j)/1.e5)**(2./7.) - 273.15
            es = 610.78*exp(17.0809*tc/(234.175+tc))
            e = rh(i,k,j)*es
            qv(i,k,j) = 0.622*e/(p_new(i,k,j)-e)
         END DO
      END DO
   END DO


END SUBROUTINE adjust_tempqv

SUBROUTINE input_terrain_rsmas ( grid ,                        &
                           ids , ide , jds , jde , kds , kde , &
                           ims , ime , jms , jme , kms , kme , &
                           ips , ipe , jps , jpe , kps , kpe )

   USE module_domain, ONLY : domain
   IMPLICIT NONE
   TYPE ( domain ) :: grid

   INTEGER , INTENT(IN)                       :: ids , ide , jds , jde , kds , kde , &
                                                 ims , ime , jms , jme , kms , kme , &
                                                 ips , ipe , jps , jpe , kps , kpe

   LOGICAL, EXTERNAL ::  wrf_dm_on_monitor

   INTEGER :: i , j , k , myproc
   INTEGER, DIMENSION(256) :: ipath  ! array for integer coded ascii for passing path down to get_terrain
   CHARACTER*256 :: message, message2
   CHARACTER*256 :: rsmas_data_path

#if DM_PARALLEL
! Local globally sized arrays
   REAL , DIMENSION(ids:ide,jds:jde) :: ht_g, xlat_g, xlon_g
#endif

   CALL wrf_get_myproc ( myproc )

#if 0
CALL domain_clock_get ( grid, current_timestr=message2 )
WRITE ( message , FMT = '(A," HT before ",I3)' ) TRIM(message2), grid%id
write(30+myproc,*)ipe-ips+1,jpe-jps+1,trim(message)
do j = jps,jpe
do i = ips,ipe
write(30+myproc,*)grid%ht(i,j)
enddo
enddo
#endif

   CALL nl_get_rsmas_data_path(1,rsmas_data_path)
   do i = 1, LEN(TRIM(rsmas_data_path))
      ipath(i) = ICHAR(rsmas_data_path(i:i))
   enddo

#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )

   CALL wrf_patch_to_global_real ( grid%xlat , xlat_g , grid%domdesc, ' ' , 'xy' ,       &
                                   ids, ide-1 , jds , jde-1 , 1 , 1 , &
                                   ims, ime   , jms , jme   , 1 , 1 , &
                                   ips, ipe   , jps , jpe   , 1 , 1   )
   CALL wrf_patch_to_global_real ( grid%xlong , xlon_g , grid%domdesc, ' ' , 'xy' ,       &
                                   ids, ide-1 , jds , jde-1 , 1 , 1 , &
                                   ims, ime   , jms , jme   , 1 , 1 , &
                                   ips, ipe   , jps , jpe   , 1 , 1   )

   IF ( wrf_dm_on_monitor() ) THEN
     CALL get_terrain ( grid%dx/1000., xlat_g(ids:ide,jds:jde), xlon_g(ids:ide,jds:jde), ht_g(ids:ide,jds:jde), &
                        ide-ids+1,jde-jds+1,ide-ids+1,jde-jds+1, ipath, LEN(TRIM(rsmas_data_path)) )
     WHERE ( ht_g(ids:ide,jds:jde) < -1000. ) ht_g(ids:ide,jds:jde) = 0.
   ENDIF

   CALL wrf_global_to_patch_real ( ht_g , grid%ht , grid%domdesc, ' ' , 'xy' ,         &
                                   ids, ide-1 , jds , jde-1 , 1 , 1 , &
                                   ims, ime   , jms , jme   , 1 , 1 , &
                                   ips, ipe   , jps , jpe   , 1 , 1   )
#else

   CALL get_terrain ( grid%dx/1000., grid%xlat(ids:ide,jds:jde), grid%xlong(ids:ide,jds:jde), grid%ht(ids:ide,jds:jde), &
                       ide-ids+1,jde-jds+1,ide-ids+1,jde-jds+1, ipath, LEN(TRIM(rsmas_data_path)) )
   WHERE ( grid%ht(ids:ide,jds:jde) < -1000. ) grid%ht(ids:ide,jds:jde) = 0.

#endif

#if 0
CALL domain_clock_get ( grid, current_timestr=message2 )
WRITE ( message , FMT = '(A," HT after ",I3)' ) TRIM(message2), grid%id
write(30+myproc,*)ipe-ips+1,jpe-jps+1,trim(message)
do j = jps,jpe
do i = ips,ipe
write(30+myproc,*)grid%ht(i,j)
enddo
enddo
#endif

END SUBROUTINE input_terrain_rsmas

SUBROUTINE update_after_feedback_em ( grid  &
!
#include "dummy_new_args.inc"
!
                 )
!
! perform core specific updates, exchanges after
! model feedback  (called from med_feedback_domain) -John
!

! Driver layer modules
   USE module_domain, ONLY : domain, get_ijk_from_grid
   USE module_configure
   USE module_driver_constants
   USE module_machine
   USE module_tiles
#ifdef DM_PARALLEL
   USE module_dm, ONLY : ntasks, ntasks_x, ntasks_y, itrace, local_communicator, mytask
   USE module_comm_dm, ONLY : HALO_EM_FEEDBACK_sub
#else
   USE module_dm
#endif
   USE module_bc
! Mediation layer modules
! Registry generated module
   USE module_state_description

   IMPLICIT NONE

   !  Subroutine interface block.

   TYPE(domain) , TARGET         :: grid

   !  Definitions of dummy arguments
#include "dummy_new_decl.inc"

   INTEGER                         :: ids , ide , jds , jde , kds , kde , &
                                      ims , ime , jms , jme , kms , kme , &
                                      ips , ipe , jps , jpe , kps , kpe

  CALL wrf_debug( 500, "entering update_after_feedback_em" )

!  Obtain dimension information stored in the grid data structure.
  CALL get_ijk_from_grid (  grid ,                   &
                            ids, ide, jds, jde, kds, kde,    &
                            ims, ime, jms, jme, kms, kme,    &
                            ips, ipe, jps, jpe, kps, kpe    )

  CALL wrf_debug( 500, "before HALO_EM_FEEDBACK.inc in update_after_feedback_em" )
#ifdef DM_PARALLEL
#include "HALO_EM_FEEDBACK.inc"
#endif
  CALL wrf_debug( 500, "leaving update_after_feedback_em" )

END SUBROUTINE update_after_feedback_em

SUBROUTINE compute_vcoord_1d_coeffs ( ht, etac, znw, &
                                      hybrid_opt, &
                                      r_d, g, p1000mb, &
                                      p_top, p00, t00, a, &
                                      ids, ide, jds, jde, kds, kde, &
                                      ims, ime, jms, jme, kms, kme, &
                                      its, ite, jts, jte, kts, kte, &
                                      znu, &
                                      c1f, c2f, c3f, c4f, &
                                      c1h, c2h, c3h, c4h )

   IMPLICIT NONE

   !  Input 

   INTEGER ,                      INTENT(IN ) :: hybrid_opt
   REAL    ,                      INTENT(IN ) :: etac
   REAL    ,                      INTENT(IN ) :: r_d, g, p1000mb, p_top, p00, t00, a
   INTEGER ,                      INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
                                                 ims, ime, jms, jme, kms, kme, &
                                                 its, ite, jts, jte, kts, kte
   REAL    , DIMENSION(kms:kme) ,           INTENT(IN ) :: znw
   REAL    , DIMENSION(ims:ime , jms:jme) , INTENT(IN ) :: ht

   !  Output

   REAL    , DIMENSION(kms:kme) , INTENT(OUT) :: znu
   REAL    , DIMENSION(kms:kme) , INTENT(OUT) :: c1f, c2f, c3f, c4f, &
                                                 c1h, c2h, c3h, c4h

   !  Local vars

   INTEGER :: i, j, k
   REAL    :: B1, B2, B3, B4, B5
   REAL    :: p_surf, press_above, press_below
   CHARACTER (LEN=256) :: a_message

   DO k=kts, kte
      IF      ( hybrid_opt .EQ. 0 ) THEN
         c3f(k) = znw(k)
      ELSE IF ( hybrid_opt .EQ. 1 ) THEN
         c3f(k) = znw(k)
      ELSE IF ( hybrid_opt .EQ. 2 ) THEN
         B1 = 2. * etac**2 * ( 1. - etac )
         B2 = -etac * ( 4. - 3. * etac - etac**3 )
         B3 = 2. * ( 1. - etac**3 )
         B4 = - ( 1. - etac**2 )
         B5 = (1.-etac)**4
         c3f(k) = ( B1 + B2*znw(k) + B3*znw(k)**2 + B4*znw(k)**3 ) / B5
         IF ( znw(k) .LT. etac ) THEN
            c3f(k) = 0.
         END IF
         IF ( k .EQ. kds ) THEN
            c3f(k) = 1.
         ELSE IF ( k .EQ. kde ) THEN
            c3f(k) = 0.
         END IF
      ELSE IF ( hybrid_opt .EQ. 3 ) THEN
         c3f(k) = znw(k)*sin(0.5*3.14159*znw(k))**2
         IF      ( k .EQ. kds ) THEN
            c3f(k) = 1.
         ELSE IF ( k .EQ. kds ) THEN
            c3f(kde) = 0.
         END IF
      ELSE
         CALL wrf_message     ( 'ERROR: --- hybrid_opt' )
         CALL wrf_message     ( 'ERROR: --- hybrid_opt=0    ==> Standard WRF terrain-following coordinate' )
         CALL wrf_message     ( 'ERROR: --- hybrid_opt=1    ==> Standard WRF terrain-following coordinate, hybrid c1, c2, c3, c4' )
         CALL wrf_message     ( 'ERROR: --- hybrid_opt=2    ==> Hybrid, Klemp polynomial' )
         CALL wrf_message     ( 'ERROR: --- hybrid_opt=3    ==> Hybrid, sin^2' )
         CALL wrf_error_fatal ( 'ERROR: --- Invalid option' )
      END IF
   END DO

   !  c4 is a function of c3 and eta.

   DO k=1, kde
      c4f(k) = ( znw(k) - c3f(k) ) * ( p1000mb - p_top )
   ENDDO

   !  Now on half levels, just add up and divide by 2 (for c3h).  Use (eta-c3)*(p00-pt) for c4 on half levels.

   DO k=1, kde-1
      znu(k) = ( znw(k+1) + znw(k) ) * 0.5
      c3h(k) = ( c3f(k+1) + c3f(k) ) * 0.5
      c4h(k) = ( znu(k) - c3h(k) ) * ( p1000mb - p_top )
   ENDDO

   !  c1 = d(B)/d(eta).  We define c1f as c1 on FULL levels.  For a vertical difference,
   !  we need to use B and eta on half levels.  The k-loop ends up referring to the
   !  full levels, neglecting the top and bottom.

   DO k=kds+1, kde-1
      c1f(k) = ( c3h(k) - c3h(k-1) ) / ( znu(k) - znu(k-1) )
   ENDDO

   !  The boundary conditions to get the coefficients:
   !  1) At k=kts: define d(B)/d(eta) = 1.  This gives us the same value of B and d(B)/d(eta)
   !     when doing the sigma-only B=eta.
   !  2) At k=kte: define d(B)/d(eta) = 0.  The curve B SMOOTHLY goes to zero, and at the very
   !     top, B continues to SMOOTHLY go to zero.  Note that for almost all cases of non B=eta,
   !     B is ALREADY=ZERO at the top, so this is a reasonable BC to assume.

   c1f(kds) = 1.
   IF ( ( hybrid_opt .EQ. 0 ) .OR. ( hybrid_opt .EQ. 1 ) ) THEN
      c1f(kde) = 1.
   ELSE
      c1f(kde) = 0.
   END IF

   !  c2 = ( 1. - c1(k) ) * (p00 - pt).  There is no vertical differencing, so we can do the
   !  full kds to kde looping.

   DO k=kds, kde
      c2f(k) = ( 1. - c1f(k) ) * ( p1000mb - p_top )
   ENDDO

   !  Now on half levels for c1 and c2.  The c1h will result from the full level c3 and full
   !  level eta differences.  The c2 value use the half level c1(k).

   DO k=1, kde-1
      c1h(k) = ( c3f(k+1) - c3f(k) ) / ( znw(k+1) - znw(k) )
      c2h(k) = ( 1. - c1h(k) ) * ( p1000mb - p_top )
   ENDDO

   !  With c3 and c4 on full levels, we can verify that we are maintaining monotonically
   !  decreasing reference pressure.
   !  P_dry_ref(k)  = c3f(k) * ( psurf - p_top ) + c4f(k) + p_top

   DO j = jts, MIN(jde-1,jte)
      DO i = its, MIN(ide-1,ite)
         p_surf = p00 * EXP ( -t00/a + ( (t00/a)**2 - 2.*g*ht(i,j)/a/r_d ) **0.5 )
         DO k=1, kde-1
            press_above = c3f(k+1)*(p_surf-p_top) + c4f(k+1) + p_top
            press_below = c3f(k  )*(p_surf-p_top) + c4f(k  ) + p_top
            IF ( press_below - press_above .LE. 0.0 ) THEN
               CALL wrf_message ( '----- ERROR: The reference pressure is not monotonically decreasing' )
               CALL wrf_message ( '             This tends to be caused by very high topography')
               WRITE (a_message,*) '(i,j) = ',i,j,', topography = ',ht(i,j),' m'
               CALL wrf_message ( TRIM(a_message) )
               WRITE (a_message,*) 'k = ',k,', reference pressure = ',press_below,' Pa'
               CALL wrf_message ( TRIM(a_message) )
               WRITE (a_message,*) 'k = ',k+1,', reference pressure = ',press_above,' Pa'
               CALL wrf_message ( TRIM(a_message) )
               WRITE (a_message,*) 'In the dynamics namelist record, reduce etac from ',etac
               CALL wrf_error_fatal ( TRIM(a_message) )
            END IF
         ENDDO
      ENDDO
   ENDDO
                                      
END SUBROUTINE compute_vcoord_1d_coeffs
