/wrfv2_fire/dyn_em/module_initialize_real.F

Large files files are truncated, but you can click here to view the full file

!REAL:MODEL_LAYER:INITIALIZATION

#ifndef VERT_UNIT
!  This MODULE holds the routines which are used to perform various initializations
!  for the individual domains, specifically for the Eulerian, mass-based coordinate.

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

MODULE module_initialize_real

   USE module_bc
   USE module_configure
   USE module_domain
   USE module_io_domain
   USE module_model_constants
   USE module_state_description
   USE module_timing
   USE module_soil_pre
   USE module_date_time
   USE module_llxy
#ifdef DM_PARALLEL
   USE module_dm
   USE module_comm_dm, ONLY : &
                           HALO_EM_INIT_1_sub   &
                          ,HALO_EM_INIT_2_sub   &
                          ,HALO_EM_INIT_3_sub   &
                          ,HALO_EM_INIT_4_sub   &
                          ,HALO_EM_INIT_5_sub   &
                          ,HALO_EM_VINTERP_UV_1_sub
#endif

   REAL , SAVE :: p_top_save
   INTEGER :: internal_time_loop

CONTAINS

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

   SUBROUTINE init_domain ( grid )

      IMPLICIT NONE

      !  Input space and data.  No gridded meteorological data has been stored, though.

!     TYPE (domain), POINTER :: grid
      TYPE (domain)          :: grid

      !  Local data.

      INTEGER :: idum1, idum2

      CALL set_scalar_indices_from_config ( head_grid%id , idum1, idum2 )

        CALL init_domain_rk( grid &
!
#include "actual_new_args.inc"
!
      )
   END SUBROUTINE init_domain

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

   SUBROUTINE init_domain_rk ( grid &
!
#include "dummy_new_args.inc"
!
   )

      USE module_optional_input
      IMPLICIT NONE

      !  Input space and data.  No gridded meteorological data has been stored, though.

!     TYPE (domain), POINTER :: grid
      TYPE (domain)          :: grid

#include "dummy_new_decl.inc"

      TYPE (grid_config_rec_type)              :: config_flags

      !  Local domain indices and counters.

      INTEGER :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat
      INTEGER :: loop , num_seaice_changes

      INTEGER :: ids, ide, jds, jde, kds, kde, &
                 ims, ime, jms, jme, kms, kme, &
                 its, ite, jts, jte, kts, kte, &
                 ips, ipe, jps, jpe, kps, kpe, &
                 i, j, k

      INTEGER :: imsx, imex, jmsx, jmex, kmsx, kmex,    &
                 ipsx, ipex, jpsx, jpex, kpsx, kpex,    &
                 imsy, imey, jmsy, jmey, kmsy, kmey,    &
                 ipsy, ipey, jpsy, jpey, kpsy, kpey

      INTEGER :: ns

      !  Local data

      INTEGER :: error
      INTEGER :: im, num_3d_m, num_3d_s
      REAL    :: p_surf, p_level
      REAL    :: cof1, cof2
      REAL    :: qvf , qvf1 , qvf2 , pd_surf
      REAL    :: p00 , t00 , a , tiso
      REAL    :: hold_znw , ptemp
      REAL    :: vap_pres_mb , sat_vap_pres_mb
      LOGICAL :: were_bad

      LOGICAL :: stretch_grid, dry_sounding, debug
      INTEGER IICOUNT

      REAL :: p_top_requested , temp
      INTEGER :: num_metgrid_levels
      REAL , DIMENSION(max_eta) :: eta_levels
      REAL :: max_dz

!      INTEGER , PARAMETER :: nl_max = 1000
!      REAL , DIMENSION(nl_max) :: grid%dn

integer::oops1,oops2

      REAL    :: zap_close_levels
      INTEGER :: force_sfc_in_vinterp
      INTEGER :: interp_type , lagrange_order , extrap_type , t_extrap_type
      LOGICAL :: lowest_lev_from_sfc , use_levels_below_ground , use_surface
      LOGICAL :: we_have_tavgsfc , we_have_tsk

      INTEGER :: lev500 , loop_count
      REAL    :: zl , zu , pl , pu , z500 , dz500 , tvsfc , dpmu
      REAL    :: pfu, pfd, phm

      LOGICAL , PARAMETER :: want_full_levels = .TRUE.
      LOGICAL , PARAMETER :: want_half_levels = .FALSE.

      CHARACTER (LEN=80) :: a_message
      REAL :: max_mf
    
      !  Excluded middle.

      LOGICAL :: any_valid_points
      INTEGER :: i_valid , j_valid

!-- Carsel and Parrish [1988]
        REAL , DIMENSION(100) :: lqmi

      REAL :: t_start , t_end

      !  Dimension information stored in grid data structure.

      CALL cpu_time(t_start)
      CALL get_ijk_from_grid (  grid ,                   &
                                ids, ide, jds, jde, kds, kde,    &
                                ims, ime, jms, jme, kms, kme,    &
                                ips, ipe, jps, jpe, kps, kpe,    &
                                imsx, imex, jmsx, jmex, kmsx, kmex,    &
                                ipsx, ipex, jpsx, jpex, kpsx, kpex,    &
                                imsy, imey, jmsy, jmey, kmsy, kmey,    &
                                ipsy, ipey, jpsy, jpey, kpsy, kpey )
      its = ips ; ite = ipe ; jts = jps ; jte = jpe ; kts = kps ; kte = kpe


      CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )

      !  Send out a quick message about the time steps based on the map scale factors.

      IF ( ( internal_time_loop .EQ. 1 ) .AND. ( grid%id .EQ. 1 ) .AND. &
           ( .NOT. config_flags%map_proj .EQ. PROJ_CASSINI ) ) THEN
         max_mf = grid%msft(its,jts)
         DO j=jts,MIN(jde-1,jte)
            DO i=its,MIN(ide-1,ite)
               max_mf = MAX ( max_mf , grid%msft(i,j) )
            END DO
         END DO
#if ( defined(DM_PARALLEL)  &&   ! defined(STUBMPI) )
         max_mf = wrf_dm_max_real ( max_mf )
#endif
         WRITE ( a_message , FMT='(A,F5.2,A)' ) 'Max map factor in domain 1 = ',max_mf, &
                                                '. Scale the dt in the model accordingly.'
         CALL wrf_message ( a_message ) 
      END IF

      !  Check to see if the boundary conditions are set properly in the namelist file.
      !  This checks for sufficiency and redundancy.

      CALL boundary_condition_check( config_flags, bdyzone, error, grid%id )

      !  Some sort of "this is the first time" initialization.  Who knows.

      grid%step_number = 0
      grid%itimestep=0

      !  Pull in the info in the namelist to compare it to the input data.

      grid%real_data_init_type = model_config_rec%real_data_init_type

      !  To define the base state, we call a USER MODIFIED routine to set the three
      !  necessary constants:  p00 (sea level pressure, Pa), t00 (sea level temperature, K),
      !  and A (temperature difference, from 1000 mb to 300 mb, K).
   
      CALL const_module_initialize ( p00 , t00 , a , tiso ) 

      !  Save these constants to write out in model output file

      grid%t00  = t00
      grid%p00  = p00
      grid%tlp  = a
      grid%tiso = tiso

      !  Are there any hold-ups to us bypassing the middle of the domain?  These
      !  holdups would be situations where we need data in the middle of the domain.
      !  FOr example, if this si the first time period, we need the full domain
      !  processed for ICs.  Also, if there is some sort of gridded FDDA turned on, or
      !  if the SST update is activated, then we can't just blow off the middle of the 
      !  domain all willy-nilly.  Other cases of these hold-ups?  Sure - what if the
      !  user wants to smooth the CG topo, we need several rows and columns available.
      !  What if the lat/lon proj is used, then we need to run a spectral filter on
      !  the topo.  Both are killers when trying to ignore data in the middle of the
      !  domain.

      !  If hold_ups = .F., then there are no hold-ups to excluding the middle
      !  domain processing.  If hold_ups = .T., then there are hold-ups, and we 
      !  must process the middle of the domain.

      hold_ups = ( internal_time_loop .EQ. 1 ) .OR. &
                 ( config_flags%grid_fdda .NE. 0 ) .OR. &
                 ( config_flags%sst_update .EQ. 1 ) .OR. &
                 ( config_flags%all_ic_times ) .OR. &
                 ( config_flags%smooth_cg_topo ) .OR. &
                 ( config_flags%map_proj .EQ. PROJ_CASSINI )

      !  There are a few checks that we need to do when the input data comes in with the middle
      !  excluded by WPS.

      IF      ( flag_excluded_middle .NE. 0 ) THEN

         !  If this time period of data from WPS has the middle excluded, it had better be OK for
         !  us to have a hole.

         IF ( hold_ups ) THEN
            WRITE ( a_message,* ) 'None of the following are allowed to be TRUE : '
            CALL wrf_message ( a_message ) 
            WRITE ( a_message,* ) ' ( internal_time_loop .EQ. 1 )               ', ( internal_time_loop .EQ. 1 )
            CALL wrf_message ( a_message ) 
            WRITE ( a_message,* ) ' ( config_flags%grid_fdda .NE. 0 )           ', ( config_flags%grid_fdda .NE. 0 )
            CALL wrf_message ( a_message ) 
            WRITE ( a_message,* ) ' ( config_flags%sst_update .EQ. 1 )          ', ( config_flags%sst_update .EQ. 1 )
            CALL wrf_message ( a_message ) 
            WRITE ( a_message,* ) ' ( config_flags%all_ic_times )               ', ( config_flags%all_ic_times )
            CALL wrf_message ( a_message ) 
            WRITE ( a_message,* ) ' ( config_flags%smooth_cg_topo )             ', ( config_flags%smooth_cg_topo )
            CALL wrf_message ( a_message ) 
            WRITE ( a_message,* ) ' ( config_flags%map_proj .EQ. PROJ_CASSINI ) ', ( config_flags%map_proj .EQ. PROJ_CASSINI )
            CALL wrf_message ( a_message ) 

            WRITE ( a_message,* ) 'Problems, we cannot have excluded middle data from WPS'
            CALL wrf_error_fatal ( a_message )
         END IF

         !  Make sure that the excluded middle data from metgrid is "wide enough".  We only have to check
         !  when the excluded middle was actually used in WPS.

         IF ( config_flags%spec_bdy_width .GT. flag_excluded_middle ) THEN
            WRITE ( a_message,* ) 'The WRF &bdy_control namelist.input spec_bdy_width = ', config_flags%spec_bdy_width
            CALL wrf_message ( a_message ) 
            WRITE ( a_message,* ) 'The WPS &metgrid namelist.wps process_only_bdy width = ',flag_excluded_middle
            CALL wrf_message ( a_message ) 
            WRITE ( a_message,* ) 'WPS process_only_bdy must be >= WRF spec_bdy_width'
            CALL wrf_error_fatal ( a_message )
         END IF
      END IF
      em_width = config_flags%spec_bdy_width

      !  We need to find if there are any valid non-excluded-middle points in this
      !  tile.  If so, then we need to hang on to a valid i,j location.

      any_valid_points = .false.
      find_valid : DO j = jts,jte
         DO i = its,ite
            IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
            any_valid_points = .true.
            i_valid = i
            j_valid = j
            EXIT find_valid
         END DO 
      END DO find_valid

      !  Replace traditional seaice field with optional seaice (AFWA source)

      IF ( flag_icefrac .EQ. 1 ) THEN
         DO j=jts,MIN(jde-1,jte)
            DO i=its,MIN(ide-1,ite)
               IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
               grid%xice(i,j) = grid%icefrac_gc(i,j)
            END DO
         END DO
      END IF

      !  Fix the snow (water equivalent depth, kg/m^2) and the snowh (physical snow
      !  depth, m) fields.

      IF      ( ( flag_snow .EQ. 0 ) .AND. ( flag_snowh .EQ. 0 ) ) THEN
         DO j=jts,MIN(jde-1,jte)
            DO i=its,MIN(ide-1,ite)
               IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
               grid%snow(i,j)  = 0.
               grid%snowh(i,j) = 0.
            END DO
         END DO

      ELSE IF ( ( flag_snow .EQ. 0 ) .AND. ( flag_snowh .EQ. 1 ) ) THEN
         DO j=jts,MIN(jde-1,jte)
            DO i=its,MIN(ide-1,ite)
!              ( m -> kg/m^2 )  & ( reduce to liquid, 5:1 ratio )
               IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
               grid%snow(i,j)  = grid%snowh(i,j) * 1000. / 5.
            END DO
         END DO

      ELSE IF ( ( flag_snow .EQ. 1 ) .AND. ( flag_snowh .EQ. 0 ) ) THEN
         DO j=jts,MIN(jde-1,jte)
            DO i=its,MIN(ide-1,ite)
!              ( kg/m^2 -> m)  & ( liquid to snow depth, 5:1 ratio )
               IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
               grid%snowh(i,j) = grid%snow(i,j) / 1000. * 5.
            END DO
         END DO

      END IF

      !  For backward compatibility, we might need to assign the map factors from
      !  what they were, to what they are.

      IF      ( ( config_flags%polar ) .AND. ( flag_mf_xy .EQ. 1 ) ) THEN
         DO j=max(jds+1,jts),min(jde-1,jte)
            DO i=its,min(ide-1,ite)
               IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
               grid%msfvx_inv(i,j) = 1./grid%msfvx(i,j)
            END DO
         END DO
         IF(jts == jds) THEN
            DO i=its,ite
               IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
               grid%msfvx(i,jts) = 0.
               grid%msfvx_inv(i,jts) = 0.
            END DO
         END IF
         IF(jte == jde) THEN
            DO i=its,ite
               IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
               grid%msfvx(i,jte) = 0.
               grid%msfvx_inv(i,jte) = 0.
            END DO
         END IF
      ELSE IF ( ( config_flags%map_proj .EQ. PROJ_CASSINI ) .AND. ( flag_mf_xy .EQ. 1 ) ) THEN
         DO j=jts,jte
            DO i=its,min(ide-1,ite)
               IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
               grid%msfvx_inv(i,j) = 1./grid%msfvx(i,j)
            END DO
         END DO
      ELSE IF ( ( .NOT. config_flags%map_proj .EQ. PROJ_CASSINI ) .AND. ( flag_mf_xy .NE. 1 ) ) THEN
         DO j=jts,jte
            DO i=its,ite
               IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
               grid%msfvx(i,j) = grid%msfv(i,j)
               grid%msfvy(i,j) = grid%msfv(i,j)
               grid%msfux(i,j) = grid%msfu(i,j)
               grid%msfuy(i,j) = grid%msfu(i,j)
               grid%msftx(i,j) = grid%msft(i,j)
               grid%msfty(i,j) = grid%msft(i,j)
            ENDDO
         ENDDO
         DO j=jts,min(jde,jte)
            DO i=its,min(ide-1,ite)
               IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
               grid%msfvx_inv(i,j) = 1./grid%msfvx(i,j)
            END DO
         END DO
      ELSE IF ( ( .NOT. config_flags%map_proj .EQ. PROJ_CASSINI ) .AND. ( flag_mf_xy .EQ. 1 ) ) THEN
         IF ( grid%msfvx(its,jts) .EQ. 0 ) THEN
            CALL wrf_error_fatal ( 'Maybe you do not have the new map factors, try re-running geogrid' )
         END IF
         DO j=jts,min(jde,jte)
            DO i=its,min(ide-1,ite)
               IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
               grid%msfvx_inv(i,j) = 1./grid%msfvx(i,j)
            END DO
         END DO
      ELSE IF ( ( config_flags%map_proj .EQ. PROJ_CASSINI ) .AND. ( flag_mf_xy .NE. 1 ) ) THEN
         CALL wrf_error_fatal ( 'Neither SI data nor older metgrid data can initialize a global domain' )
      ENDIF

      !  Check to see what available surface temperatures we have.

      IF ( flag_tavgsfc .EQ. 1 ) THEN
         we_have_tavgsfc = .TRUE.
      ELSE
         we_have_tavgsfc = .FALSE.
      END IF

      IF ( flag_tsk     .EQ. 1 ) THEN
         we_have_tsk     = .TRUE.
      ELSE
         we_have_tsk     = .FALSE.
      END IF
   
      IF ( config_flags%use_tavg_for_tsk ) THEN
         IF ( we_have_tsk .OR. we_have_tavgsfc ) THEN
           !  we are OK
         ELSE
            CALL wrf_error_fatal ( 'We either need TSK or TAVGSFC, verify these fields are coming from WPS' )
         END IF
   
         !  Since we require a skin temperature in the model, we can use the average 2-m temperature if provided.
   
         IF ( we_have_tavgsfc ) THEN
            DO j=jts,min(jde,jte)
               DO i=its,min(ide-1,ite)
                  IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
                  grid%tsk(i,j) = grid%tavgsfc(i,j)
               END DO
            END DO
         END IF
      END IF

      !  Is there any vertical interpolation to do?  The "old" data comes in on the correct
      !  vertical locations already.

      IF ( flag_metgrid .EQ. 1 ) THEN  !   <----- START OF VERTICAL INTERPOLATION PART ---->

         !  If this is data from the PINTERP program, it is emulating METGRID output.
         !  One of the caveats of this data is the way that the vertical structure is
         !  handled.  We take the k=1 level and toss it (it is disposable), and we
         !  swap in the surface data.  This is done for all of the 3d fields about
         !  which we show some interest: u, v, t, rh, ght, and p.  For u, v, and rh,
         !  we assume no interesting vertical structure, and just assign the 1000 mb
         !  data.  We directly use the 2-m temp for surface temp.  We use the surface
         !  pressure field and the topography elevation for the lowest level of
         !  pressure and height, respectively.

         IF ( flag_pinterp .EQ. 1 ) THEN

            WRITE ( a_message , * ) 'Data from P_INTERP program, filling k=1 level with artificial surface fields.'
            CALL wrf_message ( a_message )
            DO j=jts,jte
               DO i=its,ite
                  IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
                  grid%u_gc(i,1,j) = grid%u_gc(i,2,j)
                  grid%v_gc(i,1,j) = grid%v_gc(i,2,j)
                  grid%rh_gc(i,1,j) = grid%rh_gc(i,2,j)
                  grid%t_gc(i,1,j) = grid%t2(i,j)
                  grid%ght_gc(i,1,j) = grid%ht(i,j)
                  grid%p_gc(i,1,j) = grid%psfc(i,j)
               END DO
            END DO
            flag_psfc = 0

         END IF

         !  Variables that are named differently between SI and WPS.

         DO j = jts, MIN(jte,jde-1)
            DO i = its, MIN(ite,ide-1)
               IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
               grid%tsk(i,j) = grid%tsk_gc(i,j)
               grid%tmn(i,j) = grid%tmn_gc(i,j)
               grid%xlat(i,j) = grid%xlat_gc(i,j)
               grid%xlong(i,j) = grid%xlong_gc(i,j)
               grid%ht(i,j) = grid%ht_gc(i,j)
            END DO
         END DO

         !  A user could request that the most coarse grid has the
         !  topography along the outer boundary smoothed.  This smoothing
         !  is similar to the coarse/nest interface.  The outer rows and
         !  cols come from the existing large scale topo, and then the
         !  next several rows/cols are a linear ramp of the large scale
         !  model and the hi-res topo from WPS.  We only do this for the
         !  coarse grid since we are going to make the interface consistent
         !  in the model betwixt the CG and FG domains.

         IF ( ( config_flags%smooth_cg_topo ) .AND. &
              ( grid%id .EQ. 1 ) .AND. &
              ( flag_soilhgt .EQ. 1) ) THEN
            CALL blend_terrain ( grid%toposoil  , grid%ht , &
                                 ids , ide , jds , jde , 1   , 1   , &
                                 ims , ime , jms , jme , 1   , 1   , &
                                 ips , ipe , jps , jpe , 1   , 1   )

         END IF

         !  Filter the input topography if this is a polar projection.

         IF ( ( config_flags%polar ) .AND. ( grid%fft_filter_lat .GT. 90 ) ) THEN
            CALL wrf_error_fatal ( 'If the polar boundary condition is used, then fft_filter_lat must be set in namelist.input' )
         END IF

         IF ( config_flags%map_proj .EQ. PROJ_CASSINI ) THEN
#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )

            !  We stick the topo and map fac in an unused 3d array. The map scale
            !  factor and computational latitude are passed along for the ride
            !  (part of the transpose process - we only do 3d arrays) to determine
            !  "how many" values are used to compute the mean.  We want a number
            !  that is consistent with the original grid resolution.


            DO j = jts, MIN(jte,jde-1)
              DO k = kts, kte
                 DO i = its, MIN(ite,ide-1)
                    grid%t_init(i,k,j) = 1.
                 END DO
              END DO
              DO i = its, MIN(ite,ide-1)
                 grid%t_init(i,1,j) = grid%ht(i,j)
                 grid%t_init(i,2,j) = grid%msftx(i,j)
                 grid%t_init(i,3,j) = grid%clat(i,j)
              END DO
            END DO

# include "XPOSE_POLAR_FILTER_TOPO_z2x.inc"

            !  Retrieve the 2d arrays for topo, map factors, and the
            !  computational latitude.

            DO j = jpsx, MIN(jpex,jde-1)
              DO i = ipsx, MIN(ipex,ide-1)
                 grid%ht_xxx(i,j)   = grid%t_xxx(i,1,j)
                 grid%mf_xxx(i,j)   = grid%t_xxx(i,2,j)
                 grid%clat_xxx(i,j) = grid%t_xxx(i,3,j)
              END DO
            END DO

            !  Get a mean topo field that is consistent with the grid
            !  distance on each computational latitude loop.

            CALL filter_topo ( grid%ht_xxx , grid%clat_xxx , grid%mf_xxx , &
                               grid%fft_filter_lat , &
                               ids, ide, jds, jde, 1 , 1 , &
                               imsx, imex, jmsx, jmex, 1, 1, &
                               ipsx, ipex, jpsx, jpex, 1, 1 )

            !  Stick the filtered topo back into the dummy 3d array to
            !  transpose it back to "all z on a patch".

            DO j = jpsx, MIN(jpex,jde-1)
              DO i = ipsx, MIN(ipex,ide-1)
                 grid%t_xxx(i,1,j) = grid%ht_xxx(i,j)
              END DO
            END DO

# include "XPOSE_POLAR_FILTER_TOPO_x2z.inc"

            !  Get the un-transposed topo data.

            DO j = jts, MIN(jte,jde-1)
              DO i = its, MIN(ite,ide-1)
                 grid%ht(i,j) = grid%t_init(i,1,j)
              END DO
            END DO
#else
            CALL filter_topo ( grid%ht , grid%clat , grid%msftx , &
                               grid%fft_filter_lat , &
                               ids, ide, jds, jde, 1,1,           &
                               ims, ime, jms, jme, 1,1,           &
                               its, ite, jts, jte, 1,1 )
#endif
         END IF

         !  If we have any input low-res surface pressure, we store it.

         IF ( flag_psfc .EQ. 1 ) THEN
            DO j = jts, MIN(jte,jde-1)
              DO i = its, MIN(ite,ide-1)
                 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
                 grid%psfc_gc(i,j) = grid%psfc(i,j)
                 grid%p_gc(i,1,j) = grid%psfc(i,j)
              END DO
            END DO
         END IF

         !  If we have the low-resolution surface elevation, stick that in the
         !  "input" locations of the 3d height.  We still have the "hi-res" topo
         !  stuck in the grid%ht array.  The grid%landmask if test is required as some sources
         !  have ZERO elevation over water (thank you very much).

         IF ( flag_soilhgt .EQ. 1) THEN
            DO j = jts, MIN(jte,jde-1)
               DO i = its, MIN(ite,ide-1)
!                 IF ( grid%landmask(i,j) .GT. 0.5 ) THEN
                     IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
                     grid%ght_gc(i,1,j) = grid%toposoil(i,j)
                     grid%ht_gc(i,j)= grid%toposoil(i,j)
!                 END IF
               END DO
           END DO
         END IF

         !  The number of vertical levels in the input data.  There is no staggering for
         !  different variables.

         num_metgrid_levels = grid%num_metgrid_levels

         !  For UM data, swap incoming extra (theta-based) pressure with the standardly
         !  named (rho-based) pressure.

         IF ( flag_ptheta .EQ. 1 ) THEN
            DO j = jts, MIN(jte,jde-1)
               DO k = 1 , num_metgrid_levels
                  DO i = its, MIN(ite,ide-1)
                     IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
                     ptemp = grid%p_gc(i,k,j)
                     grid%p_gc(i,k,j) = grid%prho_gc(i,k,j)
                     grid%prho_gc(i,k,j) = ptemp
                  END DO
               END DO
            END DO

            !  For UM data, the "surface" and the "first hybrid" level for the theta-level data fields are the same.  
            !  Average the surface (k=1) and the second hybrid level (k=num_metgrid_levels-1) to get the first hybrid
            !  layer.  We only do this for the theta-level data: pressure, temperature, specific humidity, and
            !  geopotential height (i.e. we do not modify u, v, or the rho-based pressure).

            DO j = jts, MIN(jte,jde-1)
               DO i = its, MIN(ite,ide-1)
                  IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
                  grid%  p_gc(i,num_metgrid_levels,j) = ( grid%  p_gc(i,1,j) + grid%  p_gc(i,num_metgrid_levels-1,j) ) * 0.5
                  grid%  t_gc(i,num_metgrid_levels,j) = ( grid%  t_gc(i,1,j) + grid%  t_gc(i,num_metgrid_levels-1,j) ) * 0.5
                  grid% sh_gc(i,num_metgrid_levels,j) = ( grid% sh_gc(i,1,j) + grid% sh_gc(i,num_metgrid_levels-1,j) ) * 0.5
                  grid%ght_gc(i,num_metgrid_levels,j) = ( grid%ght_gc(i,1,j) + grid%ght_gc(i,num_metgrid_levels-1,j) ) * 0.5
               END DO
            END DO
         END IF

         IF ( any_valid_points ) THEN
         !  Check for and semi-fix missing surface fields.

         IF ( grid%p_gc(i_valid,num_metgrid_levels,j_valid) .LT. grid%p_gc(i_valid,2,j_valid) ) THEN
            k = 2
         ELSE
            k = num_metgrid_levels
         END IF

         IF ( grid%t_gc(i_valid,1,j_valid) .EQ. -1.E30 ) THEN
            DO j = jts, MIN(jte,jde-1)
               DO i = its, MIN(ite,ide-1)
                  IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
                  grid%t_gc(i,1,j) = grid%t_gc(i,k,j)
               END DO
            END DO
            config_flags%use_surface = .FALSE.
            grid%use_surface = .FALSE.
            WRITE ( a_message , * ) 'Missing surface temp, replaced with closest level, use_surface set to false.'
            CALL wrf_message ( a_message ) 
         END IF

         IF ( grid%rh_gc(i_valid,1,j_valid) .EQ. -1.E30 ) THEN
            DO j = jts, MIN(jte,jde-1)
               DO i = its, MIN(ite,ide-1)
                  IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
                  grid%rh_gc(i,1,j) = grid%rh_gc(i,k,j)
               END DO
            END DO
            config_flags%use_surface = .FALSE.
            grid%use_surface = .FALSE.
            WRITE ( a_message , * ) 'Missing surface RH, replaced with closest level, use_surface set to false.'
            CALL wrf_message ( a_message ) 
         END IF

         IF ( grid%u_gc(i_valid,1,j_valid) .EQ. -1.E30 ) THEN
            DO j = jts, MIN(jte,jde-1)
               DO i = its, ite
                  IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
                  grid%u_gc(i,1,j) = grid%u_gc(i,k,j)
               END DO
            END DO
            config_flags%use_surface = .FALSE.
            grid%use_surface = .FALSE.
            WRITE ( a_message , * ) 'Missing surface u wind, replaced with closest level, use_surface set to false.'
            CALL wrf_message ( a_message ) 
         END IF

         IF ( grid%v_gc(i_valid,1,j_valid) .EQ. -1.E30 ) THEN
            DO j = jts, jte
               DO i = its, MIN(ite,ide-1)
                  IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
                  grid%v_gc(i,1,j) = grid%v_gc(i,k,j)
               END DO
            END DO
            config_flags%use_surface = .FALSE.
            grid%use_surface = .FALSE.
            WRITE ( a_message , * ) 'Missing surface v wind, replaced with closest level, use_surface set to false.'
            CALL wrf_message ( a_message ) 
         END IF

         !  Compute the mixing ratio from the input relative humidity.

         IF ( ( flag_qv .NE. 1 ) .AND. ( flag_sh .NE. 1 ) ) THEN
            IF ( grid%p_gc(i_valid,num_metgrid_levels,j_valid) .LT. grid%p_gc(i_valid,2,j_valid) ) THEN
               k = 2
            ELSE
               k = num_metgrid_levels
            END IF

            IF      ( config_flags%rh2qv_method .eq. 1 ) THEN
               CALL rh_to_mxrat1(grid%rh_gc, grid%t_gc, grid%p_gc, grid%qv_gc ,         &
                                 config_flags%rh2qv_wrt_liquid ,                        &
                                 config_flags%qv_max_p_safe ,                           &
                                 config_flags%qv_max_flag , config_flags%qv_max_value , &
                                 config_flags%qv_min_p_safe ,                           &
                                 config_flags%qv_min_flag , config_flags%qv_min_value , &
                                 ids , ide , jds , jde , 1   , num_metgrid_levels ,     &
                                 ims , ime , jms , jme , 1   , num_metgrid_levels ,     &
                                 its , ite , jts , jte , 1   , num_metgrid_levels )
            ELSE IF ( config_flags%rh2qv_method .eq. 2 ) THEN
               CALL rh_to_mxrat2(grid%rh_gc, grid%t_gc, grid%p_gc, grid%qv_gc ,         &
                                 config_flags%rh2qv_wrt_liquid ,                        &
                                 config_flags%qv_max_p_safe ,                           &
                                 config_flags%qv_max_flag , config_flags%qv_max_value , &
                                 config_flags%qv_min_p_safe ,                           &
                                 config_flags%qv_min_flag , config_flags%qv_min_value , &
                                 ids , ide , jds , jde , 1   , num_metgrid_levels ,     &
                                 ims , ime , jms , jme , 1   , num_metgrid_levels ,     &
                                 its , ite , jts , jte , 1   , num_metgrid_levels )
            END IF


         ELSE IF ( flag_sh .EQ. 1 ) THEN
            IF ( grid%p_gc(i_valid,num_metgrid_levels,j_valid) .LT. grid%p_gc(i_valid,2,j_valid) ) THEN
               k = 2
            ELSE
               k = num_metgrid_levels
            END IF
            IF ( grid%sh_gc(i_valid,kts,j_valid) .LT. 1.e-6 ) THEN
               DO j = jts, MIN(jte,jde-1)
                  DO i = its, MIN(ite,ide-1)
                     IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
                     grid%sh_gc(i,1,j) = grid%sh_gc(i,k,j)
                  END DO
               END DO
            END IF

            DO j = jts, MIN(jte,jde-1)
               DO k = 1 , num_metgrid_levels
                  DO i = its, MIN(ite,ide-1)
                     IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
                     grid%qv_gc(i,k,j) = grid%sh_gc(i,k,j) /( 1. - grid%sh_gc(i,k,j) )
                     sat_vap_pres_mb = 0.6112*10.*EXP(17.67*(grid%t_gc(i,k,j)-273.15)/(grid%t_gc(i,k,j)-29.65))
                     vap_pres_mb = grid%qv_gc(i,k,j) * grid%p_gc(i,k,j)/100. / (grid%qv_gc(i,k,j) + 0.622 )
                     IF ( sat_vap_pres_mb .GT. 0 ) THEN
                        grid%rh_gc(i,k,j) = ( vap_pres_mb / sat_vap_pres_mb ) * 100.
                     ELSE
                        grid%rh_gc(i,k,j) = 0.
                     END IF
                  END DO
               END DO
            END DO

         END IF

         !  Some data sets do not provide a 3d geopotential height field.  

         IF ( grid%ght_gc(i_valid,grid%num_metgrid_levels/2,j_valid) .LT. 1 ) THEN
            DO j = jts, MIN(jte,jde-1)
               DO k = kts+1 , grid%num_metgrid_levels
                  DO i = its, MIN(ite,ide-1)
                     IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
                     grid%ght_gc(i,k,j) = grid%ght_gc(i,k-1,j) - &
                        R_d / g * 0.5 * ( grid%t_gc(i,k  ,j) * ( 1 + 0.608 * grid%qv_gc(i,k  ,j) ) +   &
                                          grid%t_gc(i,k-1,j) * ( 1 + 0.608 * grid%qv_gc(i,k-1,j) ) ) * &
                        LOG ( grid%p_gc(i,k,j) / grid%p_gc(i,k-1,j) )
                  END DO
               END DO
            END DO
         END IF

         !  If the pressure levels in the middle of the atmosphere are upside down, then
         !  this is hybrid data.  Computing the new surface pressure should use sfcprs2.

         IF ( grid%p_gc(i_valid,num_metgrid_levels/2,j_valid) .LT. grid%p_gc(i_valid,num_metgrid_levels/2+1,j_valid) ) THEN
            config_flags%sfcp_to_sfcp = .TRUE.
         END IF
         END IF

         !  Assign surface fields with original input values.  If this is hybrid data,
         !  the values are not exactly representative.  However - this is only for
         !  plotting purposes and such at the 0h of the forecast, so we are not all that
         !  worried.

         DO j = jts, min(jde-1,jte)
            DO i = its, min(ide,ite)
               IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
               grid%u10(i,j)=grid%u_gc(i,1,j)
            END DO
         END DO

         DO j = jts, min(jde,jte)
            DO i = its, min(ide-1,ite)
               IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
               grid%v10(i,j)=grid%v_gc(i,1,j)
            END DO
         END DO

         DO j = jts, min(jde-1,jte)
            DO i = its, min(ide-1,ite)
               IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
               grid%t2(i,j)=grid%t_gc(i,1,j)
            END DO
         END DO

         IF ( flag_qv .EQ. 1 ) THEN
            DO j = jts, min(jde-1,jte)
               DO i = its, min(ide-1,ite)
                  IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
                  grid%q2(i,j)=grid%qv_gc(i,1,j)
               END DO
            END DO
         END IF

         !  The requested ptop for real data cases.

         p_top_requested = grid%p_top_requested

         !  Compute the top pressure, grid%p_top.  For isobaric data, this is just the
         !  top level.  For the generalized vertical coordinate data, we find the
         !  max pressure on the top level.  We have to be careful of two things:
         !  1) the value has to be communicated, 2) the value can not increase
         !  at subsequent times from the initial value.

         IF ( internal_time_loop .EQ. 1 ) THEN
            CALL find_p_top ( grid%p_gc , grid%p_top , &
                              ids , ide , jds , jde , 1   , num_metgrid_levels , &
                              ims , ime , jms , jme , 1   , num_metgrid_levels , &
                              its , ite , jts , jte , 1   , num_metgrid_levels )

#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
            grid%p_top = wrf_dm_max_real ( grid%p_top )
#endif

            !  Compare the requested grid%p_top with the value available from the input data.

            IF ( p_top_requested .LT. grid%p_top ) THEN
               print *,'p_top_requested = ',p_top_requested
               print *,'allowable grid%p_top in data   = ',grid%p_top
               CALL wrf_error_fatal ( 'p_top_requested < grid%p_top possible from data' )
            END IF

            !  The grid%p_top valus is the max of what is available from the data and the
            !  requested value.  We have already compared <, so grid%p_top is directly set to
            !  the value in the namelist.

            grid%p_top = p_top_requested

            !  For subsequent times, we have to remember what the grid%p_top for the first
            !  time was.  Why?  If we have a generalized vert coordinate, the grid%p_top value
            !  could fluctuate.

            p_top_save = grid%p_top

         ELSE
            CALL find_p_top ( grid%p_gc , grid%p_top , &
                              ids , ide , jds , jde , 1   , num_metgrid_levels , &
                              ims , ime , jms , jme , 1   , num_metgrid_levels , &
                              its , ite , jts , jte , 1   , num_metgrid_levels )

#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
            grid%p_top = wrf_dm_max_real ( grid%p_top )
#endif
            IF ( grid%p_top .GT. p_top_save ) THEN
               print *,'grid%p_top from last time period = ',p_top_save
               print *,'grid%p_top from this time period = ',grid%p_top
               CALL wrf_error_fatal ( 'grid%p_top > previous value' )
            END IF
            grid%p_top = p_top_save
         ENDIF

         !  Get the monthly values interpolated to the current date for the traditional monthly
         !  fields of green-ness fraction and background albedo.

         CALL monthly_interp_to_date ( grid%greenfrac , current_date , grid%vegfra , &
                                       ids , ide , jds , jde , kds , kde , &
                                       ims , ime , jms , jme , kms , kme , &
                                       its , ite , jts , jte , kts , kte )

         CALL monthly_interp_to_date ( grid%albedo12m , current_date , grid%albbck , &
                                       ids , ide , jds , jde , kds , kde , &
                                       ims , ime , jms , jme , kms , kme , &
                                       its , ite , jts , jte , kts , kte )

         !  Get the min/max of each i,j for the monthly green-ness fraction.

         CALL monthly_min_max ( grid%greenfrac , grid%shdmin , grid%shdmax , &
                                ids , ide , jds , jde , kds , kde , &
                                ims , ime , jms , jme , kms , kme , &
                                its , ite , jts , jte , kts , kte )

         !  The model expects the green-ness values in percent, not fraction.

         DO j = jts, MIN(jte,jde-1)
           DO i = its, MIN(ite,ide-1)
              grid%vegfra(i,j) = grid%vegfra(i,j) * 100.
              grid%shdmax(i,j) = grid%shdmax(i,j) * 100.
              grid%shdmin(i,j) = grid%shdmin(i,j) * 100.
           END DO
         END DO

         !  The model expects the albedo fields as a fraction, not a percent.  Set the
         !  water values to 8%.

         DO j = jts, MIN(jte,jde-1)
            DO i = its, MIN(ite,ide-1)
               IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
               grid%albbck(i,j) = grid%albbck(i,j) / 100.
               grid%snoalb(i,j) = grid%snoalb(i,j) / 100.
               IF ( grid%landmask(i,j) .LT. 0.5 ) THEN
                  grid%albbck(i,j) = 0.08
                  grid%snoalb(i,j) = 0.08
               END IF
            END DO
         END DO

         !  Two ways to get the surface pressure.  1) If we have the low-res input surface
         !  pressure and the low-res topography, then we can do a simple hydrostatic
         !  relation.  2) Otherwise we compute the surface pressure from the sea-level
         !  pressure.
         !  Note that on output, grid%psfc is now hi-res.  The low-res surface pressure and
         !  elevation are grid%psfc_gc and grid%ht_gc (same as grid%ght_gc(k=1)).

         IF      ( ( flag_psfc    .EQ. 1 ) .AND. &
                   ( flag_soilhgt .EQ. 1 ) .AND. &
                   ( flag_slp     .EQ. 1 ) .AND. &
                   ( .NOT. config_flags%sfcp_to_sfcp ) ) THEN
            WRITE(a_message,FMT='(A)') 'Using sfcprs3 to compute psfc'
            CALL wrf_message ( a_message )
            CALL sfcprs3(grid%ght_gc, grid%p_gc, grid%ht, &
                         grid%pslv_gc, grid%psfc, &
                         ids , ide , jds , jde , 1   , num_metgrid_levels , &
                         ims , ime , jms , jme , 1   , num_metgrid_levels , &
                         its , ite , jts , jte , 1   , num_metgrid_levels )
         ELSE IF ( ( flag_psfc    .EQ. 1 ) .AND. &
                   ( flag_soilhgt .EQ. 1 ) .AND. &
                   ( config_flags%sfcp_to_sfcp ) ) THEN
            WRITE(a_message,FMT='(A)') 'Using sfcprs2 to compute psfc'
            CALL wrf_message ( a_message )
            CALL sfcprs2(grid%t_gc, grid%qv_gc, grid%ght_gc, grid%psfc_gc, grid%ht, &
                         grid%tavgsfc, grid%p_gc, grid%psfc, we_have_tavgsfc, &
                         ids , ide , jds , jde , 1   , num_metgrid_levels , &
                         ims , ime , jms , jme , 1   , num_metgrid_levels , &
                         its , ite , jts , jte , 1   , num_metgrid_levels )
         ELSE IF ( flag_slp     .EQ. 1 ) THEN
            WRITE(a_message,FMT='(A)') 'Using sfcprs  to compute psfc'
            CALL wrf_message ( a_message )
            CALL sfcprs (grid%t_gc, grid%qv_gc, grid%ght_gc, grid%pslv_gc, grid%ht, &
                         grid%tavgsfc, grid%p_gc, grid%psfc, we_have_tavgsfc, &
                         ids , ide , jds , jde , 1   , num_metgrid_levels , &
                         ims , ime , jms , jme , 1   , num_metgrid_levels , &
                         its , ite , jts , jte , 1   , num_metgrid_levels )
         ELSE
            WRITE(a_message,FMT='(3(A,I2),A,L1)') 'ERROR in psfc: flag_psfc = ',flag_psfc, &
                                               ', flag_soilhgt = ',flag_soilhgt , &
                                               ', flag_slp = ',flag_slp , & 
                                               ', sfcp_to_sfcp = ',config_flags%sfcp_to_sfcp 
            CALL wrf_message ( a_message ) 
            CALL wrf_error_fatal ( 'not enough info for a p sfc computation' )
         END IF

         !  If we have no input surface pressure, we'd better stick something in there.

         IF ( flag_psfc .NE. 1 ) THEN
            DO j = jts, MIN(jte,jde-1)
              DO i = its, MIN(ite,ide-1)
                 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
                 grid%psfc_gc(i,j) = grid%psfc(i,j)
                 grid%p_gc(i,1,j) = grid%psfc(i,j)
              END DO
            END DO
         END IF

         !  Integrate the mixing ratio to get the vapor pressure.

         CALL integ_moist ( grid%qv_gc , grid%p_gc , grid%pd_gc , grid%t_gc , grid%ght_gc , grid%intq_gc , &
                            ids , ide , jds , jde , 1   , num_metgrid_levels , &
                            ims , ime , jms , jme , 1   , num_metgrid_levels , &
                            its , ite , jts , jte , 1   , num_metgrid_levels )

         !  If this is UM data, the same moisture removed from the "theta" level pressure data can 
         !  be removed from the "rho" level pressures.  This is an approximation.  We'll revisit to
         !  see if this is a bad idea.

         IF ( flag_ptheta .EQ. 1 ) THEN
            DO j = jts, MIN(jte,jde-1)
               DO k = num_metgrid_levels-1 , 1 , -1
                  DO i = its, MIN(ite,ide-1)
                     IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
                     ptemp = ((grid%p_gc(i,k,j) - grid%pd_gc(i,k,j)) + (grid%p_gc(i,k+1,j) - grid%pd_gc(i,k+1,j)))/2
                     grid%pdrho_gc(i,k,j) = grid%prho_gc(i,k,j) - ptemp
                  END DO
               END DO
            END DO
         END IF


         !  Compute the difference between the dry, total surface pressure (input) and the
         !  dry top pressure (constant).

         CALL p_dts ( grid%mu0 , grid%intq_gc , grid%psfc , grid%p_top , &
                      ids , ide , jds , jde , 1   , num_metgrid_levels , &
                      ims , ime , jms , jme , 1   , num_metgrid_levels , &
                      its , ite , jts , jte , 1   , num_metgrid_levels )

         !  Compute the dry, hydrostatic surface pressure.

         CALL p_dhs ( grid%pdhs , grid%ht , p00 , t00 , a , &
                      ids , ide , jds , jde , kds , kde , &
                      ims , ime , jms , jme , kms , kme , &
                      its , ite , jts , jte , kts , kte )

         !  Compute the eta levels if not defined already.

         IF ( grid%znw(1) .NE. 1.0 ) THEN

            eta_levels(1:kde) = model_config_rec%eta_levels(1:kde)
            max_dz            = model_config_rec%max_dz

            CALL compute_eta ( grid%znw , &
                               eta_levels , max_eta , max_dz , &
                               grid%p_top , g , p00 , cvpm , a , r_d , cp , t00 , p1000mb , t0 , tiso , &
                               ids , ide , jds , jde , kds , kde , &
                               ims , ime , jms , jme , kms , kme , &
                               its , ite , jts , jte , kts , kte )
         END IF

         IF ( config_flags%interp_theta ) THEN

            !  The input field is temperature, we want potential temp.

            CALL t_to_theta ( grid%t_gc , grid%p_gc , p00 , &
                              ids , ide , jds , jde , 1   , num_metgrid_levels , &
                              ims , ime , jms , jme , 1   , num_metgrid_levels , &
                              its , ite , jts , jte , 1   , num_metgrid_levels )
         END IF

         IF ( flag_slp .EQ. 1 ) THEN

            !  On the eta surfaces, compute the dry pressure = mu eta, stored in
            !  grid%pb, since it is a pressure, and we don't need another kms:kme 3d
            !  array floating around.  The grid%pb array is re-computed as the base pressure
            !  later after the vertical interpolations are complete.

            CALL p_dry ( grid%mu0 , grid%znw , grid%p_top , grid%pb , want_full_levels , &
                         ids , ide , jds , jde , kds , kde , &
                         ims , ime , jms , jme , kms , kme , &
                         its , ite , jts , jte , kts , kte )

            !  All of the vertical interpolations are done in dry-pressure space.  The
            !  input data has had the moisture removed (grid%pd_gc).  The target levels (grid%pb)
            !  had the vapor pressure removed from the surface pressure, then they were
            !  scaled by the eta levels.

            interp_type = 2
            lagrange_order = grid%lagrange_order
            lowest_lev_from_sfc = .FALSE.
            use_levels_below_ground = .TRUE.
            use_surface = .TRUE.
            zap_close_levels = grid%zap_close_levels
            force_sfc_in_vinterp = 0
            t_extrap_type = grid%t_extrap_type
            extrap_type = 1

            !  For the height field, the lowest level pressure is the slp (approximately "dry").  The
            !  lowest level of the input height field (to be associated with slp) then is an array
            !  of zeros.

            DO j = jts, MIN(jte,jde-1)
               DO i = its, MIN(ite,ide-1)
                  grid%psfc_gc(i,j) = grid%pd_gc(i,1,j)
                  grid%pd_gc(i,1,j) = grid%pslv_gc(i,j) - ( grid%p_gc(i,1,j) - grid%pd_gc(i,1,j) )
                  grid%ht_gc(i,j) = grid%ght_gc(i,1,j)
                  grid%ght_gc(i,1,j) = 0.
               END DO
            END DO

            CALL vert_interp ( grid%ght_gc , grid%pd_gc , grid%ph0 , grid%pb , &
                               num_metgrid_levels , 'Z' , &
                               interp_type , lagrange_order , extrap_type , &
                               lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
                               zap_close_levels , force_sfc_in_vinterp , &
                               ids , ide , jds , jde , kds , kde , &
                               ims , ime , jms , jme , kms , kme , &
                               its , ite , jts , jte , kts , kte )

            !  Put things back to normal.

            DO j = jts, MIN(jte,jde-1)
               DO i = its, MIN(ite,ide-1)
                  IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
                  grid%pd_gc(i,1,j) = grid%psfc_gc(i,j)
                  grid%ght_gc(i,1,j) = grid%ht_gc(i,j)
               END DO
            END DO

         END IF

         !  Now the rest of the variables on half-levels to inteprolate.

         CALL p_dry ( grid%mu0 , grid%znw , grid%p_top , grid%pb , want_half_levels , &
                      ids , ide , jds , jde , kds , kde , &
                      ims , ime , jms , jme , kms , kme , &
                      its , ite , jts , jte , kts , kte )

         interp_type = grid%interp_type
         lagrange_order = grid%lagrange_order
         lowest_lev_from_sfc = grid%lowest_lev_from_sfc
         use_levels_below_ground = grid%use_levels_below_ground
         use_surface = grid%use_surface
         zap_close_levels = grid%zap_close_levels
         force_sfc_in_vinterp = grid%force_sfc_in_vinterp
         t_extrap_type = grid%t_extrap_type
         extrap_type = grid%extrap_type

         !  Interpolate RH, diagnose Qv later when have temp and pressure.  Tempo…