!
!##################################################################
!##################################################################
!######                                                      ######
!######               SUBROUTINE READNAMELIST                ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!

SUBROUTINE readnamelist(progopt,hinfmt,bdybasfn,hisfile,nhisfile,       & 2,51
     nx,ny,nz,nzsoil,nstyps,nprocx_in,nprocy_in,ncompressx,ncompressy,  &
     use_arps_grid,nx_wrf,ny_wrf,nz_wrf,zlevels_wrf,ptop,               &
     mapproj_wrf,sclfct_wrf,lattru_wrf,lontru_wrf,                      &
     ctrlat_wrf,ctrlon_wrf,dx_wrf,dy_wrf,dt,base_temp,                  &
     sfcinitopt,static_dir,start_date,silwt,wvln,                       &
     create_bdy,mgrdbas,tintv_bdyin,spec_bdy_width,                     &
     diff_opt,km_opt,khdif,kvdif,mp_physics,ra_lw_physics,              &
     ra_sw_physics,sf_sfclay_physics,sf_surface_physics,bl_pbl_physics, &
     cu_physics, nprocx_wrf,nprocy_wrf,iorder,korder,io_form,           &
     create_namelist,wrfnamelist,istatus)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  To read and propagate namelist input.
!    
!-----------------------------------------------------------------------
!
!  AUTHOR: Yunheng Wang
!  09/15/2005
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER, INTENT(IN)  :: progopt     ! 0 = ARPS2WRF
                                      ! 1 = WRFSTATIC
  INTEGER, INTENT(OUT) :: istatus

!----------------------------------------------------------------------
!
! ARPS grid variables
!
!---------------------------------------------------------------------

  INTEGER, INTENT(OUT) :: nx,ny,nz          ! Grid dimensions for ARPS.
  INTEGER, INTENT(OUT) :: nzsoil            ! Soil levels 
  INTEGER, INTENT(OUT) :: nstyps            ! Maximum number of soil types.

  INTEGER, PARAMETER :: nhisfile_max = 100
  INTEGER, PARAMETER :: max_vertical_levels = 100

  CHARACTER(LEN=*), INTENT(OUT) :: hisfile(nhisfile_max)
  CHARACTER(LEN=*), INTENT(OUT) :: bdybasfn(nhisfile_max)
  INTEGER,          INTENT(OUT) :: hinfmt
  INTEGER,          INTENT(OUT) :: nhisfile

  CHARACTER(LEN=256) :: grdbasfn
  INTEGER            :: lengbf

  ! hisfile(1) and bdybasfn(1) are for WRF input file
  !      They can be any ARPS history dumps including ADAS output, 
  !      ARPS output or EXT2ARPS output
  !
  ! hisfile(2:nhisfil), bdybasfn(2:nhisfile) are for 
  !      WRF lateral bounday files. They can be either ARPS history 
  !      dumps or EXT2ARPS outputs

  CHARACTER(LEN=4), PARAMETER :: finfmt(11) =                           &
                          (/'.bin','.asc','.hdf','.pak','.svi','.bn2',  &
                            '.net','.npk','.gad','.grb','.v5d'  /)
!
!-----------------------------------------------------------------------
!
!  ARPS include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'mp.inc'
!
!-----------------------------------------------------------------------
!
! Variables for mpi jobs
!
!-----------------------------------------------------------------------

  INTEGER, PARAMETER   :: fzone_arps = 3, fzone_wrf = 1
  INTEGER, INTENT(OUT) :: ncompressx, ncompressy ! compression in x and y direction:
                                    ! ncompressx=nprocx_in/nproc_x
                                    ! ncompressy=nprocy_in/nproc_y

  INTEGER, INTENT(OUT) :: nprocx_in, nprocy_in

!-----------------------------------------------------------------------
!
!  Namelist definitions for ARPS2WRF.input
!
!     sfcdt              Specify surface characteristics
!     bdyspc             Obtain boundary input files (ARPS format)
!     wrf_grid           Define WRF horizontal and vertical grid
!     interp_options     Choose interpolation scheme
!     wrf_opts           WRF options from namelist.input
!     output             Ouput options
!
!-----------------------------------------------------------------------
!
  INTEGER, INTENT(OUT) :: nx_wrf      ! = nx-2 if the same grid as ARPS 
  INTEGER, INTENT(OUT) :: ny_wrf      ! = ny-2 if the same grid as ARPS
  INTEGER, INTENT(OUT) :: nz_wrf      ! = nz-2 if the same grid as ARPS
                                      ! All are staggered values
                                      
  REAL,    INTENT(OUT) :: lattru_wrf(2)  ! array of true latitude of WRF map projection
  REAL,    INTENT(OUT) :: lontru_wrf     ! true longitude of WRF map projection
                                      ! = trulon_wrf
 
! Namelist variable declaration

  CHARACTER(LEN=5),   INTENT(OUT) :: sfcinitopt      ! either "ARPS" or "WRFSI"
  CHARACTER(LEN=256), INTENT(OUT) :: static_dir
  CHARACTER(LEN=19),  INTENT(OUT) :: start_date
  REAL                            :: silavwt_parm_wrf,toptwvl_parm_wrf
  REAL,               INTENT(OUT) :: silwt,wvln

  INTEGER, INTENT(OUT) :: create_bdy      ! Create WRF boundary file
  INTEGER, INTENT(OUT) :: create_namelist ! Dump WRF namelist.input
  INTEGER, INTENT(OUT) :: use_arps_grid   ! Use ARPS horizontal grid as WRF grid

  CHARACTER(LEN=80)    :: bdyfheader      ! ARPS boundary input file header
  INTEGER              :: tbgn_bdyin      ! ARPS boundary begin time (in seconds)
  INTEGER, INTENT(OUT) :: tintv_bdyin     ! ARPS boundary file interval (in seconds)
  INTEGER              :: tend_bdyin      ! Last ARPS boundary file time
  INTEGER, INTENT(OUT) :: mgrdbas         ! Options for grid base file
                                       ! = 0 share same grid base as initial state file
                                       ! = 1 All ARPS boundary files share one grd base 
                                       !     file but it is difference from the inital 
                                       !     base file as specified using grdbasfn
                                       ! = 2 Each file has its own grid base file
 
  CHARACTER(LEN=256), INTENT(OUT) :: wrfnamelist  ! file name for WRF namelist.input
 
  INTEGER, INTENT(OUT) :: mapproj_wrf     ! Type of map projection in WRF model grid
                             ! modproj = 1  Polar Stereographic
                             !              projection.
                             !         = 2  Mercator projection.
                             !         = 3  Lambert projection.
 
  REAL,    INTENT(OUT) :: sclfct_wrf      ! Map scale factor.
                             ! Distance on map, between two latitudes
                             ! trulat1 and trulat2,
                             ! is = (Distance on earth)*sclfct.
                             ! For ARPS model runs,
                             ! generally this is 1.0

  REAL              :: trulat1_wrf, trulat2_wrf, trulon_wrf
                             ! 1st, 2nd real true latitude and true longitude
                             ! of WRF map projection
  REAL, INTENT(OUT) :: ctrlat_wrf      ! Center latitude of WRF model domain (deg. N)
  REAL, INTENT(OUT) :: ctrlon_wrf      ! Center longitude of WRF model domain (deg. E)
  REAL, INTENT(OUT) :: dx_wrf          ! WRF Grid spacing in x-direction 
  REAL, INTENT(OUT) :: dy_wrf          ! WRF Grid spacing in y-direction 

  INTEGER           :: vertgrd_opt     ! WRF sigma level scheme
  REAL, INTENT(OUT) :: ptop            ! WRF atmosphere top pressure in Pascal
  REAL              :: pbot
  REAL, INTENT(OUT) :: zlevels_wrf(max_vertical_levels)
                             ! WRF mass levels from 1.0 at surfact to
                             ! 0.0 at atmosphere top

  INTEGER, INTENT(OUT) :: iorder          ! order of polynomial for horizontal 
                             ! interpolation (1 or 2)
  INTEGER, INTENT(OUT) :: korder          ! vertical interpolation order (1 or 2)

  INTEGER :: dyn_opt         ! WRF dynamics option
                             ! only works for = 2 Eulerian mass coordinate
  INTEGER, INTENT(OUT) :: diff_opt        ! WRF diffusion option
  INTEGER, INTENT(OUT) :: km_opt          ! WRF eddy coefficient option
  REAL,    INTENT(OUT) :: khdif           ! Horizontal diffusion constant (m^2/s)
  REAL,    INTENT(OUT) :: kvdif           ! Vertical diffusion constant (m^2/s)
  INTEGER, INTENT(OUT) :: mp_physics      ! WRF microphysics options
                             != 2 Lin et al. scheme 
                             !   (QVAPOR,QRAIN,QSNOW,QCLOUD,QICE,QGRAUP)
                             != 3 NCEP 3-class simple ice scheme
                             !   (QVAPOR,QCLOUD,QICE,QRAIN,QSNOW)
                             != 4 NCEP 5-class scheme
                             !   (QVAPOR,QCLOUD,QICE,QRAIN,QSNOW)
                             != 5 Ferrier (new Eta) microphysics
                             !   (QVAPOR,QCLOUD)
  INTEGER, INTENT(OUT) :: ra_lw_physics   ! Longwave radiaiton option
  INTEGER, INTENT(OUT) :: ra_sw_physics   ! Shortwave radiation option
  INTEGER, INTENT(OUT) :: sf_sfclay_physics  ! WRF surface-layer option
  INTEGER, INTENT(OUT) :: sf_surface_physics ! WRF land-surface option
                                ! = 0 no land-surface
                                !     (DO NOT use)
                                ! = 1 Thermal diffusion scheme
                                !     (nzsoil_wrf = 5)
                                ! = 2 OSU land-surface model
                                !     (nzsoil_wrf = 4)
                                ! = 3 Do not use
                                !     (nzsoil_wrf = 6)
  INTEGER, INTENT(OUT) :: bl_pbl_physics   ! boundary-layer option
  INTEGER, INTENT(OUT) :: cu_physics       ! cumulus option
  REAL,    INTENT(OUT) :: dt               ! time-step for advection
  INTEGER, INTENT(OUT) :: spec_bdy_width   ! number of rows for specified boundary values nudging
  REAL,    INTENT(OUT) :: base_temp
  INTEGER, INTENT(OUT) :: nprocx_wrf      ! Number of X direction processors for WRF run
  INTEGER, INTENT(OUT) :: nprocy_wrf      ! Number of Y direction processors for WRF run
  INTEGER, INTENT(OUT) :: io_form

  NAMELIST /message_passing/nproc_x, nproc_y, readsplit, nprocx_in, nprocy_in  
  NAMELIST /sfcdt/ sfcinitopt,ternopt,sfcfmt,sfcdtfl,                   &
                   static_dir,silavwt_parm_wrf,toptwvl_parm_wrf,start_date

  NAMELIST /bdyspc/ create_bdy,bdyfheader,tbgn_bdyin,                   &
                    tintv_bdyin,tend_bdyin,mgrdbas

  NAMELIST /wrf_grid/ use_arps_grid,nx_wrf,ny_wrf,                      &
                      mapproj_wrf, sclfct_wrf,                          &
                      trulat1_wrf, trulat2_wrf,trulon_wrf,              &
                      ctrlat_wrf, ctrlon_wrf,                           &
                      dx_wrf, dy_wrf, ptop,                             &
                      vertgrd_opt,nz_wrf,pbot,zlevels_wrf


  NAMELIST /interp_options/ iorder, korder

  NAMELIST /wrf_opts/ diff_opt, km_opt, khdif, kvdif,                   &
                      mp_physics, ra_lw_physics, ra_sw_physics,         &
                      sf_sfclay_physics, sf_surface_physics,            &
                      bl_pbl_physics, cu_physics, base_temp, dt,        &
                      spec_bdy_width, nprocx_wrf, nprocy_wrf

  NAMELIST /output/ dirname,readyfl,io_form,                            &
                    create_namelist,wrfnamelist

!-----------------------------------------------------------------------
!
!  Misc local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k,n,ifile, nlevels
  INTEGER :: lenstr, ireturn

  LOGICAL :: fexist

  REAL    :: sigma1,sigma2
  REAL    :: plevel1,plevel2
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!
! Read mpi related block
!
!-----------------------------------------------------------------------

  IF(myproc == 0) THEN
    READ(5,message_passing)
    WRITE(6,'(a)')'Namelist message_passing was successfully read.'
  END IF
  CALL mpupdatei(nproc_x,1)
  CALL mpupdatei(nproc_y,1)
  CALL mpupdatei(readsplit,1)

  IF(readsplit > 0) THEN
    nprocx_in  = nproc_x
    nprocy_in  = nproc_y
  END IF
  CALL mpupdatei(nprocx_in,1)
  CALL mpupdatei(nprocy_in,1)

  ncompressx = nprocx_in/nproc_x
  ncompressy = nprocy_in/nproc_y

  IF (mp_opt > 0 .AND. (                              &
       (MOD(nprocx_in,nproc_x) /= 0 .OR. MOD(nprocy_in,nproc_y) /= 0)   &
       .OR. (ncompressx < 1 .OR. ncompressy < 1) )  ) THEN
    IF (myproc == 0) WRITE(6,'(3x,a/,2(3x,2(a,I2)/))')                  &
      'nprocx_in (nprocy_in) must be a multiplier of nproc_x(nproc_y)', &
      'nprocx_in = ',nprocx_in, 'nprocy_in = ',nprocy_in,               &
      'nproc_x = ',  nproc_x,   'nproc_y = ', nproc_y
    CALL arpsstop('unmatched dimension size.',1);
  END IF

  IF (mp_opt == 0) THEN
    ncompressx = 1
    ncompressy = 1
    nproc_x = 1
    nproc_y = 1
    nprocx_in = 1
    nprocy_in = 1
    myproc = 0
    loc_x = 1
    loc_y = 1
    readsplit = 0
  ELSE
    CALL mpinit_var
  END IF
!
!-----------------------------------------------------------------------
!
!  Get the names of the input data files.
!
!-----------------------------------------------------------------------
!
  nhisfile  = 1

  IF(myproc == 0) THEN
    CALL get_input_file_names(hinfmt,grdbasfn,hisfile,nhisfile)

    IF(nhisfile > 1) THEN
      PRINT *, 'This program only processes one initilized data for WRF model.'
      PRINT *, 'Check the namelist file arps2wrf.input.'
      CALL arpsstop('Wrong namelist parameter.',1)
    END IF

    bdybasfn(1) = grdbasfn
    lengbf = len_trim(grdbasfn)

  END IF
  CALL mpupdatei(hinfmt,1)

  n = 1
  IF (progopt == 1) n = 0

!----------------------------------------------------------------------
!
!  Get surface characteristics file options
!
!----------------------------------------------------------------------
!
  sfcinitopt(1:5) = ' '
  sfcfmt          = 1
  sfcdtfl(1:128)  = ' '
  ternopt         = 0
  static_dir(1:128) = ' '
  start_date        = '1998-05-25_00:00:00'
  silavwt_parm_wrf  = 0.0
  toptwvl_parm_wrf  = 2.0

  IF (myproc == 0) THEN
    READ(5,sfcdt)

    IF (progopt == 0) THEN
      IF (sfcinitopt == 'ARPS' .AND. mp_opt > 0 .AND. readsplit == 0) THEN
        WRITE(grdbasfn,'(2a,2I2.2)') TRIM(sfcdtfl),'_',loc_x,loc_y
      ELSE
        WRITE(grdbasfn,'(a)') TRIM(sfcdtfl)
      END IF

      INQUIRE(FILE=trim(grdbasfn), EXIST = fexist )
      IF(.NOT. fexist) THEN
        WRITE(6,*) 'The file ',TRIM(grdbasfn),' you specified does not exist.'
        CALL arpsstop('File does not exist.',1)
      END IF
    END IF
  END IF

  lenstr = LEN_TRIM(static_dir)
  IF (static_dir(lenstr:lenstr) /= '/') THEN
    lenstr = lenstr + 1
    static_dir(lenstr:lenstr) = '/'
  END IF

  silwt = silavwt_parm_wrf
  wvln  = toptwvl_parm_wrf

  CALL mpupdatei(sfcfmt,1)
  CALL mpupdatec(sfcinitopt,5)
  CALL mpupdatec(sfcdtfl,128)

!----------------------------------------------------------------------
!
!  Get boundary file specifications
!
!----------------------------------------------------------------------
!
  create_bdy = 0
  bdyfheader = './eta25may1998'
  tbgn_bdyin = 10800
  tintv_bdyin= 10800
  tend_bdyin = 21600
  mgrdbas    = 0

  IF (myproc == 0) THEN
    READ(5,bdyspc)
    IF(create_bdy == 1 .AND. progopt == 0) THEN
      IF(tintv_bdyin < 1) THEN
        WRITE(6,'(/a,I6,a/)') 'ERROR: Boudnary interval (tintv_bdyin =',  &
                          tintv_bdyin,') is not correct. Terminating ...'
        CALL arpsstop('Wrong namelist input parameters.',1)
      END IF

      DO i = tbgn_bdyin,tend_bdyin,tintv_bdyin
        nhisfile = nhisfile + 1
        WRITE(hisfile(nhisfile),'(2a,I6.6)')                            &
                   TRIM(bdyfheader),finfmt(hinfmt),i

        IF (mp_opt > 0 .AND. readsplit <= 0 ) THEN
          WRITE(grdbasfn,'(a,a,2i2.2)') TRIM(hisfile(nhisfile)),'_',loc_x,loc_y
        ELSE
          WRITE(grdbasfn,'(a)') hisfile(nhisfile)   ! used as temporary string
        END IF
        INQUIRE(FILE=trim(grdbasfn), EXIST = fexist )
        IF(.NOT. fexist) THEN
          WRITE(6,'(/1x,3a,I6,a/)') 'WARING: The file ',TRIM(grdbasfn),  &
                  ' does not exist. Boundary file at ',i,' was skipped.'
          nhisfile = nhisfile - 1
          !CYCLE
          CALL arpsstop('File does not exist.',1)
        END IF

        IF(mgrdbas == 1 .AND. nhisfile == 2) THEN
          WRITE(bdybasfn(nhisfile),'(3a)')                              &
                   TRIM(bdyfheader),finfmt(hinfmt),'grdbas'
          n = 2      ! maximum index when checking the existance of grid and base file
        ELSE IF(mgrdbas > 1) THEN
          WRITE(bdybasfn(nhisfile),'(3a,I2.2)')                         &
                   TRIM(bdyfheader),finfmt(hinfmt),'grdbas.',nhisfile-1
          n = nhisfile
        END IF
      END DO

    END IF   ! create_bdy == 1

    !
    ! only do check for the root processor because nproc_x may 
    ! not equal nprocx_in, etc.
    !
    DO i = 1,n
      IF (mp_opt > 0 .AND. readsplit <= 0 ) THEN
        WRITE(grdbasfn,'(a,a,2i2.2)') TRIM(bdybasfn(i)),'_',loc_x,loc_y
      ELSE
        WRITE(grdbasfn,'(a)') TRIM(bdybasfn(i))   ! used as temporary string
      END IF
      INQUIRE(FILE=trim(grdbasfn), EXIST = fexist )
      IF(.NOT. fexist) THEN
        WRITE(6,*) 'The file ',TRIM(grdbasfn),' does not exist.'
        CALL arpsstop('File does not exist.',1)
      END IF
    END DO

  END IF ! myproc == 0
  CALL mpupdatei(n,1)
  CALL mpupdatei(nhisfile,1)
  CALL mpupdatec(bdybasfn,nhisfile_max*132)    ! does not include processor info.
  CALL mpupdatec(hisfile,nhisfile_max*132)

  CALL mpupdatei(create_bdy,1)
  CALL mpupdatei(mgrdbas, 1)
  CALL mpupdatei(tintv_bdyin,1)

!-----------------------------------------------------------------------
!
!  Get WRF grid options
!
!-----------------------------------------------------------------------
!
  use_arps_grid = 0
  nx_wrf  = nx-2
  ny_wrf  = ny-2

  ptop    = 5000
  pbot    = 101300
  vertgrd_opt = 0
  nz_wrf  = 31

  zlevels_wrf = -1.0

  IF (myproc == 0) THEN
    READ(5,wrf_grid)
 
    IF(use_arps_grid == 1) THEN
      WRITE(6,'(/,3(1x,a/))')   &
  '***** You have chosen to use ARPS grid in the data  as your  ******', &
  '***** WRF grid. Please note that the two fake points in each ******', &
  '***** direction of ARPS grid are not used.                   ******'  

!-----------------------------------------------------------------------
!
! Now get dimension from ARPS file
!
!-----------------------------------------------------------------------

      IF(mp_opt > 0 .AND. readsplit <= 0) THEN
        WRITE(grdbasfn,'(a,a,2i2.2)') bdybasfn(1)(1:lengbf),'_',loc_x,loc_y
        lengbf = lengbf + 5
      END IF
      CALL get_dims_from_data(hinfmt,grdbasfn(1:lengbf),                &
                              nx,ny,nz,nzsoil,nstyps,ireturn)

      IF( ireturn /= 0 ) THEN
        PRINT*,'Problem occured when trying to get dimensions from data.'
        PRINT*,'Program stopped.'
        CALL arpsstop('get_dims_from_data error.',1)
      END IF

      IF (mp_opt > 0) THEN
        IF( readsplit > 0 ) THEN
          IF( MOD(nx-fzone_arps,nproc_x) /= 0 .OR.                      &
              MOD(ny-fzone_arps,nproc_y) /= 0) THEN
            WRITE(6,'(a/,a/,4(a,i5))')      &
             'The specification of nproc_x or nproc_y is not matched with nx or ny.',&
             'nx-3 and ny-3 must be multiples of nproc_x and nproc_y respectively.', &
             'nx = ', nx, ' ny = ', ny, ' nproc_x = ',nproc_x, ' nproc_y = ',nproc_y
            nx = 0
            ny = 0
          ELSE
            nx = (nx-fzone_arps)/nproc_x + fzone_arps
            ny = (ny-fzone_arps)/nproc_y + fzone_arps
          END IF
        ELSE
          nx = (nx-fzone_arps)*ncompressx + fzone_arps
          ny = (ny-fzone_arps)*ncompressy + fzone_arps
        END IF
      END IF

      WRITE(6,'(5(a,i5))') 'nx=',nx,', ny=',ny,', nz=',nz,                &
                           ', nzsoil=',nzsoil,', nstyps=',nstyps

      nx_wrf  = nx - 2
      ny_wrf  = ny - 2
    ELSE IF (mp_opt > 0) THEN
      WRITE(6,'(1x,2a,/,a,/,2a,/,a,/,a)') 'WARING: At present, ',       &
                'arps2wrf_mpi only works when WRF horizontal grid and', &
                ' ARPS horizontal grid are the same.', 'Please set ',   &
                'use_arps_grid = 1 in arps2wrf.input.',                 &
                'Or use no-mpi version of arps2wrf.',                   &
                ' Program stopping ...'
      CALL arpsstop('MPI mode does not work.',1)
    END IF

    IF(vertgrd_opt == 0) THEN

!-----------------------------------------------------------------------
!
! Check the validation of vertical level specifications
!
!-----------------------------------------------------------------------

      nz_wrf = 0

      ! Make sure first value is 1.0
      IF (zlevels_wrf(1) /= 1.0) THEN
        WRITE(6,'(A,F6.3)') 'Bad first level: ',zlevels_wrf(1)
        WRITE(6,'(A)')      'Mass coordinate must range from 1.0 at '     &
                            //'surface to 0.0 at top.'
        CALL arpsstop('Bad WRF vertical levels.',1)
      END IF
      nlevels = 1
      ! Make sure things are decreasing to 0.0 from bottom to top.
      level_check_loop:                                                   &
      DO k = 2, max_vertical_levels
        ! See if this is a valid level, if so increment nlevels
        ! and perform additional QC checks
        IF (zlevels_wrf(k) >= 0.) THEN
          ! Check for decreasing value of META (after 2nd value found)
          IF (zlevels_wrf(k) >= zlevels_wrf(k-1)) THEN
            PRINT '(A,I2,A,I2)', 'Level ',k, 'is >= level ' , k-1
            PRINT '(A)', 'Mass must be listed in descending order!'
            CALL arpsstop('Bad WRF vertical levels.',1)
          ELSE
            nlevels = nlevels + 1
          END IF
        ELSE
          IF (nlevels < 2) THEN
            PRINT '(A)', 'Not enough levels specified'
            CALL arpsstop('Bad WRF vertical levels.',1)
          ELSE
            nz_wrf = nlevels
            ! Check to make sure top level is 0.0
            IF (zlevels_wrf(nz_wrf) /= 0.) THEN
              PRINT '(A,F6.3)', 'Bad top level value: ',zlevels_wrf(nz_wrf)
              PRINT '(A)', 'Top level must be 0.0 for Mass Eta.'
              CALL arpsstop('Bad WRF vertical levels.',1)
            END IF
            EXIT level_check_loop
          END IF
        END IF
      END DO level_check_loop

    ELSE IF (vertgrd_opt == 1) THEN
      IF(nz_wrf < 15) nz_wrf = 15
      DO k = 1,nz_wrf
        zlevels_wrf(k) = (nz_wrf-k) / (nz_wrf - 1.)
      END DO
    ELSE IF (vertgrd_opt == 2) THEN
      IF(nz_wrf < 15) nz_wrf = 15
      DO k = 1,nz_wrf,1
        zlevels_wrf(k) = SQRT( (nz_wrf-k)/(nz_wrf-1.) )
      END DO
    ELSE IF (vertgrd_opt == 3) THEN
      IF(nz_wrf < 15) nz_wrf = 15
      nlevels = 1

      DO k = 1, nz_wrf
        sigma1 =  (k-1) / (nz_wrf-1.)
        sigma2 =  SQRT( (k-1) /(nz_wrf-1.) )

        plevel1 = sigma1 * (pbot - ptop) + ptop
        plevel2 = sigma2 * (pbot - ptop) + ptop
  
        IF (plevel1 < 33333.33) THEN
          zlevels_wrf(nlevels) = sigma1
          nlevels = nlevels + 1
        END IF
        IF (plevel2 > 33333.33) THEN
          zlevels_wrf(nlevels) = sigma2
          nlevels = nlevels + 1
        END IF
      END DO  

      nz_wrf = nlevels - 1
      WRITE(6,'(a,a,I3)') 'WARNING: the actual number of vertical',     &
                        ' levels has been changed to ',nz_wrf
      !
      ! Sort the levels in decrease order (bubble sort for simplicity)
      !
      DO k = 1,nz_wrf
        DO nlevels = 1, nz_wrf-k
          IF (zlevels_wrf(nlevels+1) > zlevels_wrf(nlevels) ) THEN
            sigma1 = zlevels_wrf(nlevels)
            zlevels_wrf(nlevels) = zlevels_wrf(nlevels+1)
            zlevels_wrf(nlevels+1) = sigma1
          END IF
        END DO
      END DO

    END IF 

    !
    ! Ouput vertical levels for testing
    !
    WRITE(6,'(/2(a,I3))') 'Number of vertical levels is ',nz_wrf,       &
                       '. vertical sigma schems is ',vertgrd_opt
    WRITE(6,FMT='(1x,/,a)',ADVANCE='NO') ' levels = '
    DO nlevels = 1, nz_wrf
      IF ( nlevels /= 1 .AND. MOD(nlevels,5) == 1 ) WRITE(6,FMT='(10x)',ADVANCE='NO') 
      WRITE(6,FMT='(F8.5,A1)',ADVANCE='NO') zlevels_wrf(nlevels),', '
      IF ( MOD(nlevels,5) == 0 ) WRITE(6,*)
    END DO
    WRITE(6,*)

  END IF  ! myproc == 0

  CALL mpupdatei(nx,1)
  CALL mpupdatei(ny,1)
  CALL mpupdatei(nz,1)
  CALL mpupdatei(nzsoil,1)
  CALL mpupdatei(nstyps,1)

  IF( nx <= 0 .OR. ny <= 0 ) THEN
    CALL arpsstop('Wrong size of dimensions.',1);
  END IF

  nstyp = nstyps

  CALL mpupdatei(nx_wrf,1)
  CALL mpupdatei(ny_wrf,1)
  CALL mpupdatei(nz_wrf,1)
  CALL mpupdatei(use_arps_grid,1)
  CALL mpupdater(ptop,1)
  CALL mpupdater(zlevels_wrf,max_vertical_levels)
 
  !
  ! Map projection parameters do not need to be updated because
  ! use_arps_grid must be 1 for mpi job
  !
  lattru_wrf(1) = trulat1_wrf
  lattru_wrf(2) = trulat2_wrf
  lontru_wrf    = trulon_wrf

!-----------------------------------------------------------------------
!
!  Get WRF options (NetCDF file needs them for global attributes)
!
!-----------------------------------------------------------------------

  dyn_opt            = 2

!  IF(dyn_opt /= 2) THEN
!    WRITE(6,*) 'NOTE: ARPS2WRF only works for WRF MASS core at present'
!    WRITE(6,*) '      dyn_opt has been reseted to 2.'
!  END IF

  diff_opt           = 0 
  km_opt             = 1 
  khdif              = 0.0
  kvdif              = 0.0
  mp_physics         = 1    ! must be specified, used for moist variable determination
  ra_lw_physics      = 1
  ra_sw_physics      = 1
  sf_sfclay_physics  = 1
  sf_surface_physics = 1    ! must be specified, used for soil layer determination
  bl_pbl_physics     = 1
  cu_physics         = 1
  base_temp          = 290.
  dt                 = 40
  spec_bdy_width     = 5
  nprocx_wrf         = -1
  nprocy_wrf         = -1

  IF (myproc == 0) THEN
    READ(5,wrf_opts)

    IF(progopt == 0 .AND.                                               &
       (sf_surface_physics > 3 .OR. sf_surface_physics < 1) ) THEN
      WRITE(6,*) 'Not a valid sf_surface_physics option - ',sf_surface_physics
      WRITE(6,*) 'It must be either 1, 2 or 3.'
      CALL arpsstop('Bad WRF namelist parameters inside arps2wrf.input.',1)
    END IF
  END IF
  CALL mpupdatei(mp_physics,1)
  CALL mpupdatei(sf_surface_physics,1)
  CALL mpupdatei(spec_bdy_width,1)
  CALL mpupdater(base_temp,1)
!
!-----------------------------------------------------------------------
!
!  Get interpolation options
!
!-----------------------------------------------------------------------
!
   iorder = 3
   korder = 2

   IF (myproc == 0) READ(5,interp_options)
   CALL mpupdatei(iorder,1)
   CALL mpupdatei(korder,1)
!
!-----------------------------------------------------------------------
!
!  Get output options
!
!-----------------------------------------------------------------------
!
   dirname         = './'
   readyfl         = 0
   create_namelist = 0

   IF (myproc == 0) THEN
     READ(5,output)
     lenstr = LEN_TRIM(dirname)
     IF(lenstr > 0) THEN
       IF(dirname(lenstr:lenstr) /= '/') THEN
         dirname(lenstr+1:lenstr+1) = '/'
         lenstr = lenstr + 1
       END IF
     ELSE
       dirname = './'
     END IF
   END IF
   CALL mpupdatei(io_form,1)

!-----------------------------------------------------------------------
!
! Successfully finished with namelist input
!
!-----------------------------------------------------------------------

  istatus = 0
  
  RETURN
END SUBROUTINE readnamelist