      subroutine readinp(igrd,ngcol,ngrow,nglay,hght,phptim,hnext,
     &                   presure,ppptim,pnext,wndu,puptim,unext,wndv,
     &                   pvptim,vnext,tsrf,psptim,tsnext,temper,ptptim,
     &                   tnext,vapor,pwptim,wnext,rkvgrd,pkptim,knext,
     &                   cldwtr,ranwtr,snowtr,gplwtr,cldod,cldph,
     &                   snow,snowage,snowrat,
     &                   cigfrc,cigtim,cigwtr,cigpcp,cigent,cigdet)
      use filunit
      use grid
      use camxcom
c 
c----CAMx v7Beta6 190902
c 
c     READINP cycles through and reads all time-variant meteorological 
c     input files.  The following files are read:
c     - 3D met fields        at next update date/time, nested grids optional
c     - 2D met fields        at next update date/time, nested grids optional
c     - 3D Kv fields         at next update date/time, nested grids optional
c     - 3D Cloud/rain fields     at current date/time, nested grids optional
c     The time-rate of change for all time-interpolated fields are also
c     determined.
c 
c     Copyright 1996 - 2018
c     Ramboll
c           
c     Modifications: 
c        1/13/99   Added checks to convert input hours at 2400 to 0 on the
c                  following day
c       11/06/01   CAMx now assumes that all input file dates are in
c                  Julian format (YYJJJ) if the simulation year is 2000
c                  or greater.  Added call to DATEERR for time mismatches.
c        8/30/02   Cloud and rain files combined into one file (CAMx formatted
c                  file only), and now cloud/rain and water vapor files can be 
c                  read in for each nest
c        1/23/03   Cloud file contains new parameters:
c                  -cloud water content (cwc)
c                  -precip water content (pwc)
c                  -optical depth (cod)
c        4/2/03    Removed option for UAM-V type cloud file
c        3/3/04    Added checks for reading water files
c                  if either chemistry, dry dep or wet dep is on
c       10/12/04   Water vapor and vertical diffusivity fields are now
c                  time-interpolated
c        6/21/05   Cloud file contains new parameters:
c                  -cloud water content (cwc)
c                  -rain water content (pwr)
c                  -snow water content (pws)
c                  -graupel water content (pwg)
c                  -optical depth (cod)
c       10/24/05   Cloud water pH initialized at 5
c       12/12/05   Date/time handling improved to handle small clock drift
c        7/14/10   Added code for in-line TUV cloud adjustment
c       01/04/11   Revised for new met input format
c       02/11/11   Moved snow field from AHO to 2D met file
c       03/29/11   Support in-line TUV with aerosol optical depth
c       04/02/12   Removed RADM cloud adjustment option, cloud/aerosol 
c                  adjustments now always done with in-line TUV
c       08/08/14   Added new snow cover input fields
c       09/02/14   Added subgrid convective model
c       07/23/19   Staggered wind flag is now grid-specific
c 
c     Input arguments: 
c        igrd                grid index
c        ngcol               number of columns
c        ngrow               number of rows
c        nglay               number of layers
c        hght                current layer interface field (m)
c        presure             current pressure field (mb)
c        wndu                current u-component wind field (m/s)
c        wndv                current v-component wind field (m/s)
c        tsrf                current surface temperature field (K)
c        temper              current temperature field (K)
c        vapor               current water vapor field (ppm)
c        rkvgrd              current vertical exchange coefficient (m2/s)
c 
c     Output arguments: 
c        phptim              time-rate change of layer interface height (m/s)
c        hnext                next layer interface height field (m)
c        ppptim              time-rate change of pressure (mb/s)
c        pnext                next pressure field (mb)
c        puptim              time-rate change of u-component wind (m/s2)
c        unext                next u-component wind field (m/s)
c        pvptim              time-rate change of v-component wind (m/s2)
c        vnext                next v-component wind field (m/s)
c        psptim              time-rate change of surface temperature (K/s)
c        tsnext               next surface temperature field (m/s)
c        ptptim              time-rate change of 3-D temperature (K/s)
c        tnext                next 3-D temperature field (m/s)
c        pwptim              time-rate change of 3-D water vapor (ppm/s)
c        wnext                next 3-D water vapor field (ppm)
c        pkptim              time-rate change of 3-D vertical diff (m2/s2)
c        knext                next 3-D vertical diffusivity field (m2/s)
c        cldwtr              cloud water content (g/m3)
c        ranwtr              rain water content (g/m3)
c        snowtr              snow water content (g/m3)
c        gplwtr              graupel water content (g/m3)
c        cldod               cloud optical depth 
c        cldph               cloud water pH
c        snow                snow cover water equivalent (m)
c        snowage             snow cover age (hr)
c        snowrat             snow water accumulation rate (m/hr)
c        cigfrc              convective cloud fraction (unitless)
c        cigtim              convective cloud timescale (s)
c        cigwtr              convective cloud water (g/m3)
c        cigpcp              convective cloud precip (g/m3)
c        cigent              convective cloud entrainment (kg/m2/s)
c        cigdet              convective cloud detrainment (kg/m2/s)
c 
c     Routines Called: 
c        TIMRATES
c        CVTWIND
c        ZEROS
c        RASSGN3D
c        INTRP2D
c        INTRPV
c 
c     Called by: 
c        TSTEP_INIT
c 
      implicit none
      include 'camx.prm'
      include 'flags.inc'
c
      integer igrd,ngcol,ngrow,nglay
      real snow(ngcol,ngrow),snowage(ngcol,ngrow),snowrat(ngcol,ngrow)
      real hnext(ngcol,ngrow,nglay),pnext(ngcol,ngrow,nglay),
     &     unext(ngcol,ngrow,nglay),vnext(ngcol,ngrow,nglay),
     &     tnext(ngcol,ngrow,nglay),wnext(ngcol,ngrow,nglay),
     &     knext(ngcol,ngrow,nglay)
      real hght(ngcol,ngrow,nglay),phptim(ngcol,ngrow,nglay),
     &     presure(ngcol,ngrow,nglay),ppptim(ngcol,ngrow,nglay),
     &     wndu(ngcol,ngrow,nglay),puptim(ngcol,ngrow,nglay),
     &     wndv(ngcol,ngrow,nglay),pvptim(ngcol,ngrow,nglay),
     &     temper(ngcol,ngrow,nglay),ptptim(ngcol,ngrow,nglay),
     &     vapor(ngcol,ngrow,nglay),pwptim(ngcol,ngrow,nglay),
     &     rkvgrd(ngcol,ngrow,nglay),pkptim(ngcol,ngrow,nglay)
      real cldwtr(ngcol,ngrow,nglay),
     &     ranwtr(ngcol,ngrow,nglay),snowtr(ngcol,ngrow,nglay),
     &     gplwtr(ngcol,ngrow,nglay),cldod(ngcol,ngrow,nglay),
     &     cldph(ngcol,ngrow,nglay)
      real cigfrc(ngcol,ngrow),cigtim(ngcol,ngrow),
     &     cigwtr(ngcol,ngrow,nglay),cigpcp(ngcol,ngrow,nglay),
     &     cigent(ngcol,ngrow,nglay),cigdet(ngcol,ngrow,nglay)
      real tsrf(ngcol,ngrow),psptim(ngcol,ngrow),tsnext(ngcol,ngrow)
c
      integer hdate,iunit,idt,i,j,k,m,n,idum
      integer icigf,icigw,icigp,icige,icigd
      real whr,wmn,atim,htim,hr,snowold
      real arr3d(MXCELLS,MXCELLS,MXLAYER)
      real arr2d(MXCELLS,MXCELLS)
      character*60 string
      character*10 namevar
      character*4 namvar(10)
c
c-----Entry point
c
c-----Find the date/time at the next update interval
c
      whr = aint(time/100.)
      wmn = anint(amod(time,100.))
      atim = 100.*whr + wmn
      htim = 100.*(whr + aint((wmn + dtinp)/60.)) +
     &             amod((wmn + dtinp),60.)
      hdate = date
      if (htim.ge.2400.) then
        htim = anint(htim - 2400.)
        hdate = hdate + 1
        if( MOD(hdate,1000) .GT. 365 ) then
           if( MOD(INT(hdate/1000),4) .EQ. 0 ) then
              if( MOD(hdate,1000) .EQ. 367 )
     &           hdate = (INT(hdate/1000)+1)*1000 + 1
           else
              hdate = (INT(hdate/1000)+1)*1000 + 1
           endif
        endif
      endif
c
c-----Read 3D met file
c
      iunit = i3dmet(igrd)
      if( iunit .GT. 0 ) then
         write(string,'(A,I5)') 'Reading 3D met file for grid:',igrd
 101     continue
         read(iunit,end=7000) idt,hr
         hr = 100*aint(hr) + anint(60.*amod(hr,1.))
         if (hr.ge.2400.) then
           hr = hr - 2400.
           idt = idt + 1
           if( MOD(idt,1000) .GT. 365 ) then
              if( MOD(INT(idt/1000),4) .EQ. 0 ) then
                 if( MOD(idt,1000) .EQ. 367 )
     &              idt = (INT(idt/1000)+1)*1000 + 1
              else
                 idt = (INT(idt/1000)+1)*1000 + 1
              endif
           endif
         endif
         write(iout,'(a40,f7.0,i8.5,a,i3)')
     &        'Read 3D met file at ',hr,idt,' grid',igrd
         call flush(iout)
c
         if (idt.lt.hdate .or. (idt.eq.hdate .and. hr.lt.htim)) then
           do n = 1,n3dmet(igrd)
             do k = 1,nglay
               read(iunit)
             enddo
           enddo
           goto 101
         endif
         if (idt.gt.hdate .or. (idt.eq.hdate .and. hr.gt.htim))
     &     goto 7001
c
         do n = 1,n3dmet(igrd)
           do k = 1,nglay
            read(iunit,end=7000) idum,(namvar(m),m=1,10),
     &                           ((arr3d(i,j,k),i=1,ngcol),j=1,ngrow)
           enddo
           write(namevar,'(10a1)') (namvar(m),m=1,10)
           if (namevar.eq.'ZGRID_M') then
             do k = 1,nglay
               do j = 1,ngrow
                 do i = 1,ngcol
                   hnext(i,j,k) = arr3d(i,j,k)
                 enddo
               enddo
             enddo
           elseif (namevar.eq.'PRESS_MB') then
             do k = 1,nglay
               do j = 1,ngrow
                 do i = 1,ngcol
                   pnext(i,j,k) = arr3d(i,j,k)
                 enddo
               enddo
             enddo
           elseif (namevar.eq.'TEMP_K') then
             do k = 1,nglay
               do j = 1,ngrow
                 do i = 1,ngcol
                   tnext(i,j,k) = arr3d(i,j,k)
                 enddo
               enddo
             enddo
           elseif (namevar.eq.'HUMID_PPM') then
             do k = 1,nglay
               do j = 1,ngrow
                 do i = 1,ngcol
                   wnext(i,j,k) = arr3d(i,j,k)
                 enddo
               enddo
             enddo
           elseif (namevar.eq.'UWIND_MpS') then
             do k = 1,nglay
               do j = 1,ngrow
                 do i = 1,ngcol
                   unext(i,j,k) = arr3d(i,j,k)
                 enddo
               enddo
             enddo
           elseif (namevar.eq.'VWIND_MpS') then
             do k = 1,nglay
               do j = 1,ngrow
                 do i = 1,ngcol
                   vnext(i,j,k) = arr3d(i,j,k)
                 enddo
               enddo
             enddo
           endif
         enddo
c
c-----Convert wind vectors if they are not staggered
c
         if (.not.lstagw(igrd)) call cvtwind(ngcol,ngrow,nglay,unext,
     &                                       vnext)
c
c-----Calculate time rates of change for 3D met fields
c
         call timrates(ngcol,ngrow,nglay,hght,hnext,phptim)
         call timrates(ngcol,ngrow,nglay,presure,pnext,ppptim)
         call timrates(ngcol,ngrow,nglay,temper,tnext,ptptim)
         call timrates(ngcol,ngrow,nglay,vapor,wnext,pwptim)
         call timrates(ngcol,ngrow,nglay,wndu,unext,puptim)
         call timrates(ngcol,ngrow,nglay,wndv,vnext,pvptim)
      endif
c
c-----Read 2D met file
c
      iunit = i2dmet(igrd)
      if( iunit .GT. 0 ) then
         write(string,'(A,I5)') 'Reading 2D met file for grid:',igrd
 102     continue
         read(iunit,end=7000) idt,hr
         hr = 100*aint(hr) + anint(60.*amod(hr,1.))
         if (hr.ge.2400.) then
           hr = hr - 2400.
           idt = idt + 1
           if( MOD(idt,1000) .GT. 365 ) then
              if( MOD(INT(idt/1000),4) .EQ. 0 ) then
                 if( MOD(idt,1000) .EQ. 367 )
     &              idt = (INT(idt/1000)+1)*1000 + 1
              else
                 idt = (INT(idt/1000)+1)*1000 + 1
              endif
           endif
         endif
         write(iout,'(a40,f7.0,i8.5,a,i3)')
     &        'Read 2D met file at ',hr,idt,' grid',igrd
         call flush(iout)
c
         if (idt.lt.hdate .or. (idt.eq.hdate .and. hr.lt.htim)) then
           do n = 1,n2dmet(igrd)
             read(iunit)
           enddo
           goto 102
         endif
         if (idt.gt.hdate .or. (idt.eq.hdate .and. hr.gt.htim))
     &     goto 7001
c
         do n = 1,n2dmet(igrd)
           read(iunit,end=7000) idum,(namvar(m),m=1,10),
     &                          ((arr2d(i,j),i=1,ngcol),j=1,ngrow)
           write(namevar,'(10a1)') (namvar(m),m=1,10)
           if (namevar.eq.'TSURF_K') then
             do j = 1,ngrow
               do i = 1,ngcol
                 tsnext(i,j) = arr2d(i,j)
               enddo
             enddo
           elseif (namevar.eq.'SNOWEW_M') then
             do j = 1,ngrow
               do i = 1,ngcol
                 snowold = snow(i,j)
                 snow(i,j) = arr2d(i,j)
                 snowrat(i,j) = (snow(i,j) - snowold)/(dtinp/60.)
               enddo
             enddo
           elseif (namevar.eq.'SNOWAGE_HR') then
             do j = 1,ngrow
               do i = 1,ngcol
                 snowage(i,j) = arr2d(i,j)
               enddo
             enddo
           elseif (namevar.eq.'SNOWCVR') then
             do j = 1,ngrow
               do i = 1,ngcol
                 if (arr2d(i,j).gt.0.) then
                   snow(i,j) = .025          ! Make some assumptions
                   snowage(i,j) = 5.*24.     ! for backward compatibility
                   snowrat(i,j) = 0.
                 endif
               enddo
             enddo
           endif
         enddo
c
c-----Calculate time rates of change for 2D met fields
c
         call timrates(ngcol,ngrow,1,tsrf,tsnext,psptim)
      endif
c
c-----Read 3D VDiff file
c
      iunit = ikv(igrd)
      if( iunit .GT. 0 ) then
         write(string,'(A,I5)') 'Reading 3D VDiff file for grid:',igrd
 103     continue
         read(iunit,end=7000) idt,hr
         hr = 100*aint(hr) + anint(60.*amod(hr,1.))
         if (hr.ge.2400.) then
           hr = hr - 2400.
           idt = idt + 1
           if( MOD(idt,1000) .GT. 365 ) then
              if( MOD(INT(idt/1000),4) .EQ. 0 ) then
                 if( MOD(idt,1000) .EQ. 367 )
     &              idt = (INT(idt/1000)+1)*1000 + 1
              else
                 idt = (INT(idt/1000)+1)*1000 + 1
              endif
           endif
         endif
         write(iout,'(a40,f7.0,i8.5,a,i3)')
     &        'Read 3D VDiff file at ',hr,idt,' grid',igrd
         call flush(iout)
c
         if (idt.lt.hdate .or. (idt.eq.hdate .and. hr.lt.htim)) then
           do n = 1,nkvmet(igrd)
             do k = 1,nglay
               read(iunit)
             enddo
           enddo
           goto 103
         endif
         if (idt.gt.hdate .or. (idt.eq.hdate .and. hr.gt.htim))
     &     goto 7001
c
         do n = 1,nkvmet(igrd)
           do k = 1,nglay
            read(iunit,end=7000) idum,(namvar(m),m=1,10),
     &                           ((arr3d(i,j,k),i=1,ngcol),j=1,ngrow)
           enddo
           write(namevar,'(10a1)') (namvar(m),m=1,10)
           if (namevar.eq.'KV_M2pS') then
             do k = 1,nglay
               do j = 1,ngrow
                 do i = 1,ngcol
                   knext(i,j,k) = arr3d(i,j,k)
                 enddo
               enddo
             enddo
           endif
         enddo
c
c-----Calculate time rates of change for 3D vertical diffusivity
c
         call timrates(ngcol,ngrow,nglay,rkvgrd,knext,pkptim)
      endif
c
c-----Read 3D cloud file
c
      iunit = icld(igrd)
      if (igrd.eq.1 .and. iunit.eq.0) then
         call zeros(cldwtr,ngcol*ngrow*nglay)
         call zeros(ranwtr,ngcol*ngrow*nglay)
         call zeros(snowtr,ngcol*ngrow*nglay)
         call zeros(gplwtr,ngcol*ngrow*nglay)
         call zeros(cldod, ngcol*ngrow*nglay)
      elseif (iunit.gt.0) then
         write(string,'(A,I5)') 'Reading 3D Cloud file for grid:',igrd
 104     continue
         read(iunit,end=7000) idt,hr
         hr = 100*aint(hr) + anint(60.*amod(hr,1.))
         if (hr.ge.2400.) then
           hr = hr - 2400.
           idt = idt + 1
           if( MOD(idt,1000) .GT. 365 ) then
              if( MOD(INT(idt/1000),4) .EQ. 0 ) then
                 if( MOD(idt,1000) .EQ. 367 )
     &              idt = (INT(idt/1000)+1)*1000 + 1
              else
                 idt = (INT(idt/1000)+1)*1000 + 1
              endif
           endif
         endif
         write(iout,'(a40,f7.0,i8.5,a,i3)')
     &        'Read 3D Cloud file at ',hr,idt,' grid',igrd
         call flush(iout)
c
         if (idt.lt.date .or. (idt.eq.date .and. hr.lt.atim)) then
           do n = 1,ncldmet(igrd)
             do k = 1,nglay
               read(iunit)
             enddo
           enddo
           goto 104
         endif
         if (idt.gt.date .or. (idt.eq.date .and. hr.gt.atim))
     &     goto 7002
c
         icigf = 0 
         icigw = 0 
         icigp = 0 
         icige = 0 
         icigd = 0 
         do n = 1,ncldmet(igrd)
           do k = 1,nglay
            read(iunit,end=7000) idum,(namvar(m),m=1,10),
     &                           ((arr3d(i,j,k),i=1,ngcol),j=1,ngrow)
           enddo
           write(namevar,'(10a1)') (namvar(m),m=1,10)
           if (namevar.eq.'CLODW_GpM3') then
             do k = 1,nglay
               do j = 1,ngrow
                 do i = 1,ngcol
                   cldwtr(i,j,k) = arr3d(i,j,k)
                 enddo
               enddo
             enddo
           endif
           if (namevar.eq.'RAINW_GpM3') then
             do k = 1,nglay
               do j = 1,ngrow
                 do i = 1,ngcol
                   ranwtr(i,j,k) = arr3d(i,j,k)
                 enddo
               enddo
             enddo
           endif
           if (namevar.eq.'SNOWW_GpM3') then
             do k = 1,nglay
               do j = 1,ngrow
                 do i = 1,ngcol
                   snowtr(i,j,k) = arr3d(i,j,k)
                 enddo
               enddo
             enddo
           endif
           if (namevar.eq.'GRPLW_GpM3') then
             do k = 1,nglay
               do j = 1,ngrow
                 do i = 1,ngcol
                   gplwtr(i,j,k) = arr3d(i,j,k)
                 enddo
               enddo
             enddo
           endif
           if (namevar.eq.'CLOUDOD') then
             do k = 1,nglay
               do j = 1,ngrow
                 do i = 1,ngcol
                   cldod(i,j,k) = arr3d(i,j,k)
                 enddo
               enddo
             enddo
           endif
           if (namevar.eq.'KFFRAC' .and. lcig(igrd)) then
             icigf = 1
             do j = 1,ngrow
               do i = 1,ngcol
                 cigfrc(i,j) = arr3d(i,j,1)
               enddo
             enddo
             do j = 1,ngrow
               do i = 1,ngcol
                 cigtim(i,j) = arr3d(i,j,2)
               enddo
             enddo
           endif
           if (namevar.eq.'KFWTR_GpM3' .and. lcig(igrd)) then
             icigw = 1
             do k = 1,nglay
               do j = 1,ngrow
                 do i = 1,ngcol
                   cigwtr(i,j,k) = arr3d(i,j,k)
                 enddo
               enddo
             enddo
           endif
           if (namevar.eq.'KFPCP_GpM3' .and. lcig(igrd)) then
             icigp = 1
             do k = 1,nglay
               do j = 1,ngrow
                 do i = 1,ngcol
                   cigpcp(i,j,k) = arr3d(i,j,k)
                 enddo
               enddo
             enddo
           endif
           if (namevar.eq.'KFE_KGpM2S' .and. lcig(igrd)) then
             icige = 1
             do k = 1,nglay
               do j = 1,ngrow
                 do i = 1,ngcol
                   cigent(i,j,k) = arr3d(i,j,k)
                 enddo
               enddo
             enddo
           endif
           if (namevar.eq.'KFD_KGpM2S' .and. lcig(igrd)) then
             icigd = 1
             do k = 1,nglay
               do j = 1,ngrow
                 do i = 1,ngcol
                   cigdet(i,j,k) = arr3d(i,j,k)
                 enddo
               enddo
             enddo
           endif
         enddo
         if (lcig(igrd) .and. icigf+icigw+icigp+icige+icigd.ne.5) then
           lcig(igrd) = .false.
           write(iout,'(//,a)')'WARNING in READINP:'
           write(iout,'(a)') string
           write(iout,'(2a)')'Subgrid convective cloud data were not',
     &                       ' found on file'
           write(iout,'(2a)')'Convective treatment will not be',
     &                       ' applied for this grid'
           write(iout,*)
         endif
c
c-----Initialize cloud pH
c
         do k = 1,nglay
           do j = 1,ngrow
             do i = 1,ngcol
               cldph(i,j,k) = 5.
             enddo
           enddo
         enddo
c
      endif
c
      return
c
 7000 continue
      write(iout,'(//,a)')'ERROR in READINP:'
      write(iout,'(a)') string
      write(iout,*)'End of input file reached.  Make sure the file '
      write(iout,*)'is for the correct day and contains all hours.'
      write(iout,*)
      call camxerr()
c
 7001 continue
      write(iout,'(//,a)') 'ERROR in READINP:'
      write(iout,'(a)') string
      write(iout,'(a,f10.1,i10.5)') 
     &      'Past expected time/date ',htim,hdate
      write(iout,*)
      call camxerr()
c
 7002 continue
      write(iout,'(//,a)') 'ERROR in READINP:'
      write(iout,'(a)') string
      write(iout,'(a,f10.1,i10.5)') 
     &      'Past expected time/date ',atim,date
      write(iout,*)
      call camxerr()
c
      end
