#include <misc.h>
#include <params.h>


module inidat 1,3
!BOP
!
! !MODULE: inidat --- dynamics-physics coupling module
!
! !USES:
   use shr_kind_mod, only: r8 => shr_kind_r8
   use comsrf
   use abortutils, only: endrun

! !PUBLIC MEMBER FUNCTIONS:
   public read_inidat, copy_inidat

! !PUBLIC DATA MEMBERS:
   real(r8), allocatable :: ps_tmp(:,:)
   real(r8), allocatable :: phis_tmp(:,:)
   real(r8), allocatable :: landfrac_tmp(:,:)                  
   real(r8), allocatable :: tsocn_tmp(:,:)                  
   real(r8), allocatable :: icefrac_tmp(:,:)                  
   real(r8), allocatable :: landm_tmp(:,:)
   real(r8), allocatable :: pblht_tmp(:,:)
   real(r8), allocatable :: tpert_tmp(:,:)
   real(r8), allocatable :: qpert_tmp(:,:)
   real(r8), allocatable :: sgh_tmp(:,:)
   real(r8), allocatable :: tsice_tmp(:,:)
   real(r8), allocatable :: tsice_rad_tmp(:,:)                   
   real(r8), allocatable :: tbot_tmp(:,:)                  
   real(r8), allocatable :: tssub_tmp(:,:,:)
   real(r8), allocatable :: sicthk_tmp(:,:)
   real(r8), allocatable :: snowhice_tmp(:,:)                

   real(r8) zgsint_tmp

   logical read_tsicerad
   logical read_tbot
   logical read_pblh
   logical read_tpert
   logical read_qpert
   logical read_cloud
   logical read_qcwat
   logical read_tcwat
   logical read_lcwat

!
! !DESCRIPTION:
!
!      This module provides 
!
!      \begin{tabular}{|l|l|} \hline \hline
!        read\_inidat    &   \\ \hline
!        copy\_inidat    &   \\ \hline 
!                                \hline
!      \end{tabular}
!
! !REVISION HISTORY:
!   YY.MM.DD   ?????      Creation
!   00.06.01   Grant      First attempt at modifying for LRDC
!   01.10.01   Lin        Various revisions
!   01.01.15   Sawyer     Bug fixes for SPMD mode
!   01.03.26   Sawyer     Added ProTeX documentation
!   02.04.04   Sawyer     Removed comspe
!   03.08.31   Mirin      Eliminated many 3D temporary global arrays
!
!EOP
!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------
!BOP
! !IROUTINE: read_inidat --- read initial dataset
!
! !INTERFACE: 

   subroutine read_inidat 1,152

! !USES:
      use pmgrid
      use pspect
      use rgrid
      use commap
      use prognostics
      use physconst,    only: gravit
      use history,      only: fillvalue
      use constituents, only: pcnst, pnats, cnst_name, qmin, cnst_read_iv
      use chemistry,    only: chem_implements_cnst, chem_init_cnst
      ! TBH:  combine modules aerosol_intr and aerosols?  
      use aerosol_intr, only: aerosol_implements_cnst, aerosol_init_cnst
      use tracers, only: tracers_implements_cnst, tracers_init_cnst
      use cldcond,      only: cldcond_implements_cnst, cldcond_init_cnst
      use phys_buffer,  only: pbuf, pbuf_times, pbuf_get_fld_idx
      use phys_grid
#if ( defined SPMD )
      use mpishorthand
      use spmd_dyn, only : comm_y, comm_z
      use parutilitiesmodule, only: parcollective2d, BCSTOP
#endif

      implicit none

      include 'netcdf.inc'

!------------------------------Commons----------------------------------

#include <comctl.h>
#include <comqfl.h>
#include <comlun.h>
#include <perturb.h>

! !DESCRIPTION:
!
!   Read initial dataset and spectrally truncate as appropriate.
!
! !REVISION HISTORY:
!
!   00.06.01   Grant      First attempt at modifying for LRDC
!   00.10.01   Lin        Various revisions
!   01.01.15   Sawyer     Bug fixes for SPMD mode
!   01.03.09   Eaton      Modifications
!   01.03.26   Sawyer     Added ProTeX documentation
!   03.08.31   Mirin      Eliminated many 3D temporary global arrays
!
!EOP
!-----------------------------------------------------------------------
!BOC
!
! !LOCAL VARIABLES:

      integer i,j,k,m,lat,lon   ! grid and constituent indices
      integer ihem              ! hemisphere index
      real(r8) pdelb(plond,plev)! pressure diff between interfaces
      real(r8) pertval          ! perturbation value
      real(r8) zgssum           ! partial sums of phis
      integer ii, ic, n, mq, ml1, ml2
!
! Netcdf related variables
!
      integer lonsiz, latsiz, levsiz ! Dimension sizes
      integer londimid, levdimid, latdimid ! Dimension ID's
      integer tid                 ! Variable ID's
      integer tracid(pcnst+pnats) ! Variable ID's
      integer phisid, sghid, psid ! Variable ID's
      integer landmid
      integer pblhtid
      integer tpertid
      integer qpertid
      integer cldid
      integer qcwatid
      integer tcwatid
      integer lcwatid
#if ( ! defined COUP_CSM )
      integer ts1id, ts2id, ts3id, ts4id,tsiceid,tsice_rad_id ! Variable ID's
      integer tbotid             ! Variable ID's
#endif
#if ( defined COUP_SOM )
      integer sicid
      integer icefracid
      integer tsocnid
#endif
      integer snowhiceid           ! Variable ID's
      integer landfracid        ! Variable ID's
      integer usid, vsid
      integer strt2d(3)         ! start lon, lat, time indices for netcdf 2-d
      integer strt3d(4)         ! start lon, lev, lat, time for netcdf 3-d
      data strt2d/3*1/          ! Only index 2 will ever change
      data strt3d/4*1/          ! Only indices 2,3 will ever change

      integer cnt2d(3)          ! lon, lat, time counts for netcdf 2-d
      integer cnt3d(4)          ! lon, lat, lev, time counts for netcdf 2-d
      data cnt2d/plon,1,1/      ! 2-d arrs: Always grab only a "plon" slice
      data cnt3d/plon,plev,plat,1/ ! 3-d arrs: Always grab a full time slice

      integer ndims2d           ! number of dimensions
      integer dims2d(NF_MAX_VAR_DIMS) ! variable shape
      integer ndims3d           ! number of dimensions
      integer dims3d(NF_MAX_VAR_DIMS) ! variable shape
      integer tmptype
      integer natt, ret, attlen ! netcdf return values
      logical phis_hires        ! true => PHIS came from hi res topo
      real(r8), allocatable :: arrxyz(:,:,:)
      real(r8), allocatable :: arrxzy(:,:,:)
      real(r8), allocatable :: uv_local(:,:,:)
      real(r8), pointer, dimension(:,:,:,:) :: cld
      real(r8), pointer, dimension(:,:,:,:) :: tcwat
      real(r8), pointer, dimension(:,:,:,:) :: qcwat
      real(r8), pointer, dimension(:,:,:,:) :: lcwat

      character*(NF_MAX_NAME) tmpname
      character*256 text
      character*80 trunits      ! tracer untis

      integer istat
      integer slatid, slatdimid, slatsiz
      integer slonid, slondimid, slonsiz
      integer cnt3dus(4)        ! index counts for netcdf U staggered grid
      integer cnt3dvs(4)        ! index counts for netcdf V staggered grid

!      data cnt3dus/plon,plev,splat,1/ ! 3-d arrs: Always grab a full time slice
! SJL
      integer platm1
      parameter (platm1=plat-1)
      data cnt3dus/plon,plev,platm1,1/ ! temporary patch
      data cnt3dvs/plon,plev,plat,1/ ! 3-d arrs: Always grab a full time slice

!
!-----------------------------------------------------------------------
!     August 2003 revision described below (Mirin)
!-----------------------------------------------------------------------
! This routine has the master process reading in global arrays from disk
!   and storing them in temporary arrays. The data is then scattered to
!   the other processes in copy_inidat. Originally, all the data was
!   read in, and subsequently all of it was scattered.
! Because of the large volume of temporary global data, particularly at
!   high resolution, the procedure was modified to use fewer temporary
!   global arrays. This required interleaving reading and scattering of
!   data. Scattering of the 3D arrays (and a few others) is now part of
!   read_inidat; the remaining quantities are scattered in copy_inidat.
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
! Allocate memory for temporary arrays
!-----------------------------------------------------------------------
!
! Note if not masterproc still might need to allocate array for spmd case
! since each processor calls MPI_scatter 
!
      allocate ( ps_tmp(plond,plat), stat=istat )
      if ( istat /= 0 ) call endrun
      allocate ( arrxyz(plond,plat,plev), stat=istat )
      if ( istat /= 0 ) call endrun
      allocate ( arrxzy(plond,plev,plat), stat=istat )             
      if ( istat /= 0 ) call endrun
      allocate( uv_local(plond,beglat:endlat,beglev:endlev), stat=istat )
      if ( istat /= 0 ) call endrun
      allocate ( phis_tmp(plond,plat), stat=istat )        
      if ( istat /= 0 ) call endrun
      allocate ( landm_tmp(plond,plat), stat=istat )                  
      if ( istat /= 0 ) call endrun
      allocate ( sgh_tmp(plond,plat), stat=istat )                 
      if ( istat /= 0 ) call endrun
      allocate ( tsice_tmp(plond,plat), stat=istat )                   
      if ( istat /= 0 ) call endrun
      allocate ( tsice_rad_tmp(plond,plat), stat=istat )                   
      if ( istat /= 0 ) call endrun
      allocate ( tbot_tmp    (plond,plat), stat=istat )                
      if ( istat /= 0 ) call endrun
      allocate ( tssub_tmp(plond,plevmx,plat), stat=istat )         
      if ( istat /= 0 ) call endrun
      allocate ( sicthk_tmp(plond,plat), stat=istat )               
      if ( istat /= 0 ) call endrun
      allocate ( snowhice_tmp(plond,plat), stat=istat )                
      if ( istat /= 0 ) call endrun
      allocate ( landfrac_tmp(plond,plat), stat=istat )                
      if ( istat /= 0 ) call endrun
      allocate ( pblht_tmp(plond,plat), stat=istat )
      if ( istat /= 0 ) call endrun
      allocate ( tpert_tmp(plond,plat), stat=istat )
      if ( istat /= 0 ) call endrun
      allocate ( qpert_tmp(plond,plat), stat=istat )
      if ( istat /= 0 ) call endrun
#if ( defined COUP_SOM )
      allocate ( icefrac_tmp(plond,plat) )                
      allocate ( tsocn_tmp(plond,plat) )                
#endif
!
!-----------------------------------------------------------------------
! Read in input variables
!-----------------------------------------------------------------------

!
! logical flags to track which "extra" fields are indeed in IC file
!
      read_tsicerad = .false.
      read_tbot     = .false.
      read_pblh     = .false.
      read_tpert    = .false.
      read_qpert    = .false.
      read_cloud    = .false.
      read_qcwat    = .false.
      read_tcwat    = .false.
      read_lcwat    = .false.
      qpert_tmp(:plon,:) = 0.

      if (masterproc) then
!
! Get dimension IDs and lengths 
!
         call wrap_inq_dimid  (ncid_ini, 'lat', latdimid)
         call wrap_inq_dimlen (ncid_ini, latdimid, latsiz)
         call wrap_inq_dimid  (ncid_ini, 'lev', levdimid)
         call wrap_inq_dimlen (ncid_ini, levdimid, levsiz)
         call wrap_inq_dimid  (ncid_ini, 'lon', londimid)
         call wrap_inq_dimlen (ncid_ini, londimid, lonsiz)
         call wrap_inq_dimid  (ncid_ini, 'slat', slatdimid)
         call wrap_inq_dimlen (ncid_ini, slatdimid, slatsiz)
         call wrap_inq_dimid  (ncid_ini, 'slon', slondimid)
         call wrap_inq_dimlen (ncid_ini, slondimid, slonsiz)

!
! Get variable id's 
! Check that all tracer units are in mass mixing ratios
!
!         call wrap_inq_varid (ncid_ini, 'U'   , uid)
!         call wrap_inq_varid (ncid_ini, 'V'   , vid)

         call wrap_inq_varid (ncid_ini, 'slat', slatid)
         call wrap_inq_varid (ncid_ini, 'slon', slonid)
         call wrap_inq_varid (ncid_ini, 'US'  , usid)
         call wrap_inq_varid (ncid_ini, 'VS'  , vsid)

         call wrap_inq_varid (ncid_ini, 'T'   , tid)
         call wrap_inq_varid (ncid_ini, 'PS'  , psid)
         call wrap_inq_varid (ncid_ini, 'PHIS', phisid)
         call wrap_inq_varid (ncid_ini, 'SGH' , sghid)
         if (nf_inq_varid (ncid_ini, 'LANDM_COSLAT', landmid) /= nf_noerr) then
            write(6,*)'INIDAT: LANDM_COSLAT not found on initial dataset.'
            write(6,*)'        Need to run definesurf to create it.'
            write(6,*)'        This field became a requirement as of cam2_0_2_dev43'
            call endrun ()
         end if

         if ( nf_inq_varid(ncid_ini, 'PBLH', pblhtid) == NF_NOERR ) then
            call wrap_inq_varid (ncid_ini, 'PBLH' , pblhtid)
            read_pblh  = .true.
         end if
         if ( nf_inq_varid(ncid_ini, 'TPERT', tpertid) == NF_NOERR ) then
            call wrap_inq_varid (ncid_ini, 'TPERT', tpertid)
            read_tpert = .true.
         end if
         if ( nf_inq_varid(ncid_ini, 'QPERT', qpertid) == NF_NOERR ) then
            call wrap_inq_varid (ncid_ini, 'QPERT', qpertid)
            read_qpert = .true.
         end if
         if ( nf_inq_varid(ncid_ini, 'CLOUD', cldid) == NF_NOERR ) then
            call wrap_inq_varid (ncid_ini, 'CLOUD', cldid  )
            read_cloud = .true.
         end if
         if ( nf_inq_varid(ncid_ini, 'QCWAT', qcwatid) == NF_NOERR ) then
            call wrap_inq_varid (ncid_ini, 'QCWAT', qcwatid)
            read_qcwat = .true.
         end if
         if ( nf_inq_varid(ncid_ini, 'TCWAT', tcwatid) == NF_NOERR ) then
            call wrap_inq_varid (ncid_ini, 'TCWAT', tcwatid)
            read_tcwat = .true.
         end if
         if ( nf_inq_varid(ncid_ini, 'LCWAT', lcwatid) == NF_NOERR ) then
            call wrap_inq_varid (ncid_ini, 'LCWAT', lcwatid)
            read_lcwat = .true.
         end if

#if ( ! defined COUP_CSM )
!
! For land-fraction check if the variable name LANDFRAC is on the dataset if not assume FLAND
!
         if ( nf_inq_varid(ncid_ini, 'LANDFRAC', landfracid ) == NF_NOERR ) then
            call wrap_inq_varid (ncid_ini, 'LANDFRAC', landfracid)
         else
            call wrap_inq_varid (ncid_ini, 'FLAND', landfracid)
         end if
         if ( nf_inq_varid(ncid_ini, 'TBOT', tbotid) == NF_NOERR ) then
            call wrap_inq_varid (ncid_ini, 'TBOT'    , tbotid)
            read_tbot     = .true.
         end if
         if ( nf_inq_varid(ncid_ini, 'TSICERAD', tsice_rad_id) == NF_NOERR ) then
            call wrap_inq_varid (ncid_ini, 'TSICERAD', tsice_rad_id)
            read_tsicerad = .true.
         end if
         call wrap_inq_varid (ncid_ini, 'TSICE', tsiceid)
         call wrap_inq_varid (ncid_ini, 'TS1', ts1id)
         call wrap_inq_varid (ncid_ini, 'TS2', ts2id)
         call wrap_inq_varid (ncid_ini, 'TS3', ts3id)
         call wrap_inq_varid (ncid_ini, 'TS4', ts4id)
         call wrap_inq_varid (ncid_ini, 'SNOWHICE', snowhiceid)
#if ( defined COUP_SOM )
         call wrap_inq_varid (ncid_ini, 'SICTHK', sicid)
         call wrap_inq_varid (ncid_ini, 'ICEFRAC', icefracid)
         call wrap_inq_varid (ncid_ini, 'TSOCN', tsocnid)
#endif
#endif
!
! Guard:  Check that "Q" is on IC file
!
         if (cnst_read_iv(1) .and. &
              nf_inq_varid(ncid_ini, cnst_name(1), tracid(1) ) /= NF_NOERR ) then
            call endrun ('READ_INIDAT'//cnst_name(1)//' not found on IC file')
         end if

         do m=1,pcnst+pnats
            if (cnst_read_iv(m) .and. &
               nf_inq_varid(ncid_ini, cnst_name(m), tracid(m) ) == NF_NOERR ) then
               call wrap_inq_varid (NCID_INI,cnst_name(m), tracid(m))
               call wrap_get_att_text (NCID_INI,tracid(m),'units', trunits)
               if (trunits(1:5) .ne. 'KG/KG' .and. trunits(1:5) .ne. 'kg/kg') then
                  call endrun ('INIDAT: tracer units for tracer = '//cnst_name(m)//' must be in KG/KG')
               endif
            end if
         end do
!
! Check dimension ordering for one 2-d and one 3-d field.
! Assume other arrays of like rank will have dimensions ordered the same.
!
         call wrap_inq_var (ncid_ini, psid, tmpname, tmptype, &
                            ndims2d, dims2d, natt)
         if (dims2d(1).ne.londimid .or. dims2d(2).ne.latdimid .or. &
             ndims2d.gt.3) then
            call endrun ('INIDAT: Bad number of dims or ordering on 2d fld')
         end if
!
! Check for presence of 'from_hires' attribute to decide whether to filter
!
         ret = nf_inq_attlen (ncid_ini, phisid, 'from_hires', attlen)
         if (ret.eq.NF_NOERR .and. attlen.gt.256) then
            call endrun ('INIDAT: from_hires attribute length is too long')
         end if
         ret = nf_get_att_text (ncid_ini, phisid, 'from_hires', text)
         if (ret.eq.NF_NOERR .and. text(1:4).eq.'true') then
            phis_hires = .true.
!            write(6,*)'INIDAT: Will filter input PHIS: attribute ', &
!                      'from_hires is true'
         else
            phis_hires = .false.
!            write(6,*)'INIDAT: Will not filter input PHIS: attribute ', &
!                      'from_hires is either false or not present'
         end if
!
! Read in 2d fields.  
! For stand alone run: get surface temp and 4 (sub)surface temp fields
! For stand alone run with slab-ocean: get sea-ice thickness and snow cover
!
         do j=1,plat
            strt2d(2) = j
            if (ideal_phys .or. aqua_planet) then
               do i=1,nlon(j)
                  phis_tmp(i,j) = 0.
                  sgh_tmp (i,j) = 0.
               end do
            else
               call wrap_get_vara_realx (ncid_ini, phisid, strt2d, cnt2d, &
                                         phis_tmp(1,j))
               call wrap_get_vara_realx (ncid_ini, sghid , strt2d, cnt2d, &
                                         sgh_tmp(1,j))
            endif
            call wrap_get_vara_realx (ncid_ini, landmid, strt2d, cnt2d, &
                                      landm_tmp(1,j))
            call wrap_get_vara_realx (ncid_ini, psid, strt2d, cnt2d, &
                                      ps_tmp(1,j))
#if ( ! defined COUP_CSM )
            if (aqua_planet) then
               do i=1,nlon(j)
                  landfrac_tmp(i,j) = 0.
                  tbot_tmp(i,j) = 0.
               end do
            else
               call wrap_get_vara_realx (ncid_ini, landfracid, strt2d, cnt2d, &
                                         landfrac_tmp(1,j))
               if (read_tbot) then
                  call wrap_get_vara_realx (ncid_ini, tbotid    , strt2d, cnt2d, tbot_tmp    (1,j))
               end if
            endif
            call wrap_get_vara_realx (ncid_ini, tsiceid, strt2d, cnt2d, &
                                           tsice_tmp(1,j))
            call wrap_get_vara_realx (ncid_ini, ts1id, strt2d, cnt2d, &
                                           tssub_tmp(1,1,j))
            call wrap_get_vara_realx (ncid_ini, ts2id, strt2d, cnt2d, &
                                           tssub_tmp(1,2,j))
            call wrap_get_vara_realx (ncid_ini, ts3id, strt2d, cnt2d, &
                                           tssub_tmp(1,3,j))
            call wrap_get_vara_realx (ncid_ini, ts4id, strt2d, cnt2d, &
                                           tssub_tmp(1,4,j))
            if (read_tsicerad) then
               call wrap_get_vara_realx (ncid_ini, tsice_rad_id, strt2d, cnt2d, tsice_rad_tmp(1,j))
            end if
            if (read_pblh) then
               call wrap_get_vara_realx (ncid_ini, pblhtid, strt2d, cnt2d, pblht_tmp(1,j  ))
            end if
            if (read_tpert) then
               call wrap_get_vara_realx (ncid_ini, tpertid, strt2d, cnt2d, tpert_tmp(1,j  ))
            end if
            if (read_qpert) then
               call wrap_get_vara_realx (ncid_ini, qpertid, strt2d, cnt2d, qpert_tmp(1,j))
            end if
!
! Set sea-ice thickness and snow cover:
!
#if ( defined COUP_SOM )
            call wrap_get_vara_realx(ncid_ini, sicid, strt2d, cnt2d, sicthk_tmp(1,j))
            call wrap_get_vara_realx(ncid_ini, icefracid, strt2d, cnt2d, icefrac_tmp(1,j))
            call wrap_get_vara_realx(ncid_ini, tsocnid, strt2d, cnt2d, tsocn_tmp(1,j))
#endif
            call wrap_get_vara_realx(ncid_ini, snowhiceid, strt2d, cnt2d, snowhice_tmp(1,j))
#endif
         end do
!
! Read in 3d fields.  

! Staggered grid variables and transpose

! Read in u3

         call wrap_get_vara_realx(ncid_ini,usid, strt3d, cnt3dus, arrxzy)

         do k = 1, plev
!
! SJL: initialize j=1 because later on arrxyz will be copied to u3s using f90 array syntax
               do i = 1, plon
                  arrxyz(i,1,k) = fillvalue
               enddo
!
            do j = 1, plat-1
               do i = 1, plon
                  arrxyz(i,j+1,k) = arrxzy(i,k,j)
               enddo
            enddo
         enddo

      endif                     ! end of if-masterproc

! Scatter u3

#if ( defined SPMD )
      call scatter( arrxyz, strip3dxyz, uv_local, mpicom )
!$omp parallel do private(i,j,k)
      do k=beglev,endlev
         do j=beglat,endlat
            do i=1,plon
               u3s(i,j,k) = uv_local(i,j,k)
            enddo
         enddo
      enddo
#else
!$omp parallel do private(i, j, k)
      do k=beglev,endlev
         do j=beglat,endlat
            do i=1,plon
               u3s(i,j,k) = arrxyz(i,j,k)
            enddo
         enddo
      enddo
#endif

      if (masterproc) then

! Read in v3

         call wrap_get_vara_realx(ncid_ini,vsid, strt3d, cnt3dvs, arrxzy)

         do k = 1, plev
            do j = 1, plat
               do i = 1, plon
                  arrxyz(i,j,k) = arrxzy(i,k,j)
               enddo
            enddo
         enddo

      endif                     ! end of if-masterproc

! Scatter v3

#if ( defined SPMD )
      call scatter( arrxyz, strip3dxyz, uv_local, mpicom )
!$omp parallel do private(i,j,k)
      do k=beglev,endlev
         do j=beglat,endlat
            do i=1,plon
               v3s(i,j,k) = uv_local(i,j,k)
            enddo
         enddo
      enddo
#else
!$omp parallel do private(i, j, k)
      do k=beglev,endlev
         do j=beglat,endlat
            do i=1,plon
               v3s(i,j,k) = arrxyz(i,j,k)
            enddo
         enddo
      enddo
#endif

      if (masterproc) then

! Read in t3

         call wrap_get_vara_realx(ncid_ini, tid, strt3d, cnt3d, arrxzy)

!         
! Add random perturbation to temperature if required
!
         if (pertlim.ne.0.0) then
            write(6,*)'INIDAT: Adding random perturbation bounded by +/-', &
                      pertlim,' to initial temperature field'
            do lat=1,plat
               do k=1,plev
                  do i=1,nlon(lat)
                     call random_number (pertval)
                     pertval = 2.*pertlim*(0.5 - pertval)
                     arrxzy(i,k,lat) = arrxzy(i,k,lat)*(1. + pertval)
                  end do
               end do
            end do
         endif
!$omp parallel do private(k)
         do k = 1, plev
            call xpavg(arrxzy(1,k,   1), plon)
            call xpavg(arrxzy(1,k,plat), plon)
         enddo

      endif                     ! end of if-masterproc

! Scatter t3

#if ( defined SPMD )
      call scatter( arrxzy, strip3dxzy, t3, mpicom )
#else
!$omp parallel do private(i, j, k, ic)
      do j=beglat,endlat
         do k=beglev,endlev
            do i=1,plon
               t3(i,k,j) = arrxzy(i,k,j)
            enddo
         enddo
      enddo
#endif

! IMPORTANT - the following block of coding must be kept adjacent to the t3
!    coding, as it assumes that ARRXZY contains t3 (for .not. read_tbot)

      if (masterproc) then
         if(.not. read_tbot) then
            tbot_tmp   (:plon,:) = arrxzy(:plon,plev,:)
         endif
      endif

! Read tcwat

! IMPORTANT - the following block of coding must be kept adjacent to the t3
!    coding, as it assumes that ARRXZY contains t3 (for .not. read_tcwat)

      if (masterproc) then

         if (read_tcwat) then
            call wrap_get_vara_realx(ncid_ini, tcwatid, strt3d, cnt3d, arrxzy)
         else
            write(6,*) 'Warning:  TCWAT not found on IC file; initialized with T'
         end if

      endif

! Scatter tcwat

      m = pbuf_get_fld_idx('TCWAT')
      tcwat => pbuf(m)%fld_ptr(1,1:pcols,1:pver,begchunk:endchunk,1:pbuf_times)
      call scatter_field_to_chunk(1,plev,1,plond,arrxzy,tcwat(:,:,:,1))

! Read cld

      if (masterproc) then

         if (read_cloud) then
            call wrap_get_vara_realx(ncid_ini, cldid  , strt3d, cnt3d, arrxzy)
         else
            arrxzy    (:plon,:,:) = 0.
            write(6,*) 'Warning:  CLOUD not found on IC file; initialized to 0.'
         end if

      endif

! Scatter cld

      m = pbuf_get_fld_idx('CLD')
      cld => pbuf(m)%fld_ptr(1,1:pcols,1:pver,begchunk:endchunk,1:pbuf_times)
      call scatter_field_to_chunk(1,plev,1,plond,arrxzy,cld(:,:,:,1))

! Read qcwat

      if (masterproc) then

         if (read_qcwat) then
            call wrap_get_vara_realx(ncid_ini, qcwatid, strt3d, cnt3d, arrxzy)
         else
! Search among constituents
            mq = 0
            do m = 1, pcnst+pnats
              if (cnst_name(m) == 'Q') mq = m
            enddo
            if (mq .eq. 0) then
               call endrun ('READ_INIDAT:Q not found - exiting')
            endif
            m = mq
            write(6,*) 'Warning:  QCWAT not found on IC file; initialized with ',cnst_name(m)
            if (cnst_read_iv(m) .and. &
               nf_inq_varid(ncid_ini, cnst_name(m), tracid(m) ) == NF_NOERR ) then
               call wrap_get_vara_realx(ncid_ini, tracid(m), strt3d, cnt3d, arrxzy)
            else
! Initialize instead of input
               write(6,*) 'Warning:  Not reading ',cnst_name(m), ' from IC file.'
               arrxzy = 0.
               if (cldcond_implements_cnst(cnst_name(m))) then
                  call cldcond_init_cnst(cnst_name(m), arrxzy)
                  write(6,*) '          ', cnst_name(m), ' initialized by "cldcond_init_cnst"'
               else if (chem_implements_cnst(cnst_name(m))) then
                  call chem_init_cnst(cnst_name(m), arrxzy)
                  write(6,*) '          ', cnst_name(m), ' initialized by "chem_init_cnst"'
               else if (tracers_implements_cnst(cnst_name(m))) then
                  call tracers_init_cnst(cnst_name(m), arrxzy)
                  write(6,*) '          ', cnst_name(m), ' initialized by "tracers_init_cnst"'
               else
                  write(6,*) '          ', cnst_name(m), ' set to 0.'
               end if
            endif
            do lat=1,plat
               call qneg3('INIDAT  ', lat, nlon(lat), plon, plev, 1, &
                          qmin(m), arrxzy(:plon,:plev,lat))
            end do
            do k = 1, plev
               call xpavg(arrxzy(:,k,   1), plon)
               call xpavg(arrxzy(:,k,plat), plon)
            enddo
         end if

      endif

! Scatter qcwat

      m = pbuf_get_fld_idx('QCWAT')
      qcwat => pbuf(m)%fld_ptr(1,1:pcols,1:pver,begchunk:endchunk,1:pbuf_times)
      call scatter_field_to_chunk(1,plev,1,plond,arrxzy,qcwat(:,:,:,1))

! Read lcwat

      if (masterproc) then
         if (read_lcwat) then
            call wrap_get_vara_realx(ncid_ini, lcwatid, strt3d, cnt3d, arrxzy)
         else
! Search among constituents
            ml1 = 0
            ml2 = 0
            do m = 1, pcnst+pnats
              if (cnst_name(m) == 'CLDLIQ') ml1 = m
              if (cnst_name(m) == 'CLDICE') ml2 = m
            enddo
            if (ml1 .eq. 0) then
               call endrun ('READ_INIDAT:CLDLIQ not found - exiting')
            endif
            if (ml2 .eq. 0) then
               call endrun ('READ_INIDAT:CLDICE not found - exiting')
            endif
            m = ml1
            write(6,*) 'Warning:  LCWAT not found on IC file; ',cnst_name(m),' added to LCWAT instead'
            if (cnst_read_iv(m) .and. &
               nf_inq_varid(ncid_ini, cnst_name(m), tracid(m) ) == NF_NOERR ) then
               call wrap_get_vara_realx(ncid_ini, tracid(m), strt3d, cnt3d, arrxzy)
            else
! Initialize instead of input
               write(6,*) 'Warning:  Not reading ',cnst_name(m), ' from IC file.'
               arrxzy = 0.
               if (cldcond_implements_cnst(cnst_name(m))) then
                  call cldcond_init_cnst(cnst_name(m), arrxzy)
                  write(6,*) '          ', cnst_name(m), ' initialized by "cldcond_init_cnst"'
               else if (chem_implements_cnst(cnst_name(m))) then
                  call chem_init_cnst(cnst_name(m), arrxzy)
                  write(6,*) '          ', cnst_name(m), ' initialized by "chem_init_cnst"'
               else if (tracers_implements_cnst(cnst_name(m))) then
                  call tracers_init_cnst(cnst_name(m), arrxzy)
                  write(6,*) '          ', cnst_name(m), ' initialized by "tracers_init_cnst"'
               else
                  write(6,*) '          ', cnst_name(m), ' set to 0.'
               end if
            endif
            do lat=1,plat
               call qneg3('INIDAT  ', lat, nlon(lat), plon, plev, 1, &
                          qmin(m), arrxzy(:plon,:plev,lat))
            end do
            do k = 1, plev
               call xpavg(arrxzy(:,k,   1), plon)
               call xpavg(arrxzy(:,k,plat), plon)
            enddo

! Store CLDLIQ in arrxyz
            do k = 1, plev
              do lat = 1, plat
                do lon = 1, plon
                  arrxyz(lon,lat,k) = arrxzy(lon,k,lat)
                enddo
              enddo
            enddo

            m = ml2
            write(6,*) 'Warning:  LCWAT not found on IC file; ',cnst_name(m),' added to LCWAT instead'
            if (cnst_read_iv(m) .and. &
               nf_inq_varid(ncid_ini, cnst_name(m), tracid(m) ) == NF_NOERR ) then
               call wrap_get_vara_realx(ncid_ini, tracid(m), strt3d, cnt3d, arrxzy)
            else
! Initialize instead of input
               write(6,*) 'Warning:  Not reading ',cnst_name(m), ' from IC file.'
               arrxzy = 0.
               if (cldcond_implements_cnst(cnst_name(m))) then
                  call cldcond_init_cnst(cnst_name(m), arrxzy)
                  write(6,*) '          ', cnst_name(m), ' initialized by "cldcond_init_cnst"'
               else if (chem_implements_cnst(cnst_name(m))) then
                  call chem_init_cnst(cnst_name(m), arrxzy)
                  write(6,*) '          ', cnst_name(m), ' initialized by "chem_init_cnst"'
               else if (tracers_implements_cnst(cnst_name(m))) then
                  call tracers_init_cnst(cnst_name(m), arrxzy)
                  write(6,*) '          ', cnst_name(m), ' initialized by "tracers_init_cnst"'
               else
                  write(6,*) '          ', cnst_name(m), ' set to 0.'
               end if
            endif
            do lat=1,plat
               call qneg3('INIDAT  ', lat, nlon(lat), plon, plev, 1, &
                          qmin(m), arrxzy(:plon,:plev,lat))
            end do
            do k = 1, plev
               call xpavg(arrxzy(:,k,   1), plon)
               call xpavg(arrxzy(:,k,plat), plon)
            enddo

! Add CLDLIQ and CLDICE
            do k = 1, plev
              do lat = 1, plat
                do lon = 1, plon
                  arrxzy(lon,k,lat) = arrxyz(lon,lat,k) + arrxzy(lon,k,lat)
                enddo
              enddo
            enddo
         end if

      endif

! Scatter lcwat

      m = pbuf_get_fld_idx('LCWAT')
      lcwat => pbuf(m)%fld_ptr(1,1:pcols,1:pver,begchunk:endchunk,1:pbuf_times)
      call scatter_field_to_chunk(1,plev,1,plond,arrxzy,lcwat(:,:,:,1))

      if (pbuf_times > 1) then
         do n = 2, pbuf_times
            cld  (:,:,:,n) = cld  (:,:,:,1)
            tcwat(:,:,:,n) = tcwat(:,:,:,1)
            qcwat(:,:,:,n) = qcwat(:,:,:,1)
            lcwat(:,:,:,n) = lcwat(:,:,:,1)
         end do
      end if

      if (masterproc) then
!
! Compute integrals of mass, moisture, and geopotential height
!
!gg  Integrals of mass and moisture should be unnecessary in Lin-Rood dynamics
!gg  because they are conserved. What's left is the global geopotential...
!gg  Dunno if that's necessary or not, so I left it in.

         zgsint_tmp = 0.
         do lat = 1, plat
!              
! Accumulate average mass of atmosphere
!
            zgssum = 0.
            do i=1,nlon(lat)
               zgssum = zgssum + phis_tmp(i,lat)
            end do
            zgsint_tmp = zgsint_tmp + w(lat)*zgssum/nlon(lat)
         end do                  ! end of latitude loop
!
! Normalize average height
!
         zgsint_tmp = zgsint_tmp*.5/gravit
!
! Globally avgd sfc. partial pressure of dry air (i.e. global dry mass):
!
! SJL:
         tmass0 = 98222./gravit
!
! WS: Average pole information moved here: treat the global arrays
!
!-----------------------------------------------------------
! Average T, PS, PHIS and Q at the poles.       The initial
! conditions *should* already have these variables averaged,
! but do it here for safety -- no harm if it's already done.
!-----------------------------------------------------------

         call xpavg(phis_tmp(1, 1), plon)
         call xpavg(phis_tmp(1,plat), plon)
         call xpavg(ps_tmp(1,   1), plon)
         call xpavg(ps_tmp(1,plat), plon)

         if (ideal_phys) tmass0 = 100000./gravit
!        write(6,800) tmassf_tmp,tmass0,qmassf_tmp
!        write(6,810) zgsint_tmp
800      format('INIDAT: MASS OF INITIAL DATA BEFORE CORRECTION = ' &
                ,1p,e20.10,/,' DRY MASS WILL BE HELD = ',e20.10,/, &
                ' MASS OF MOISTURE AFTER REMOVAL OF NEGATIVES = ',e20.10) 
810      format(/69('*')/'INIDAT: Globally averaged geopotential ', &
                'height = ',f16.10,' meters'/69('*')/)

      endif                     ! end of if-masterproc

! Scatter ps and phis

#if ( defined SPMD )
      if (myid_z .eq. 0) then
         call scatter( ps_tmp, strip2d, ps, comm_y )
         call scatter( phis_tmp, strip2d, phis, comm_y )
      endif
      if (twod_decomp .eq. 1) then
         call parcollective2d( comm_z, BCSTOP, plon, endlat-beglat+1, ps ) 
         call parcollective2d( comm_z, BCSTOP, plon, endlat-beglat+1, phis ) 
      endif
#else
      ps(:,:) = ps_tmp(:,:)
      phis(:,:) = phis_tmp(:,:)
#endif

! Initialize constituents

      do m = 1, pcnst+pnats

         if (masterproc) then

            if (cnst_read_iv(m) .and. &
               nf_inq_varid(ncid_ini, cnst_name(m), tracid(m) ) == NF_NOERR ) then
               call wrap_get_vara_realx(ncid_ini, tracid(m), strt3d, cnt3d, arrxzy)
            else
               write(6,*) 'Warning:  Not reading ',cnst_name(m), ' from IC file.'
               arrxzy = 0.
               if (cldcond_implements_cnst(cnst_name(m))) then
                  call cldcond_init_cnst(cnst_name(m), arrxzy)
                  write(6,*) '          ', cnst_name(m), ' initialized by "cldcond_init_cnst"'
               else if (chem_implements_cnst(cnst_name(m))) then
                  call chem_init_cnst(cnst_name(m), arrxzy)
                  write(6,*) '          ', cnst_name(m), ' initialized by "chem_init_cnst"'
               else if (tracers_implements_cnst(cnst_name(m))) then
                  call tracers_init_cnst(cnst_name(m), arrxzy)
                  write(6,*) '          ', cnst_name(m), ' initialized by "tracers_init_cnst"'
               else if (aerosol_implements_cnst(cnst_name(m))) then
                  call aerosol_init_cnst(cnst_name(m), arrxzy)
                  write(6,*) '          ', cnst_name(m), ' initialized by "aerosol_init_cnst"'
               else
                  write(6,*) '          ', cnst_name(m), ' set to 0.'
               end if
            endif

!$omp parallel do private(lat)
            do lat=1,plat
               call qneg3('INIDAT  ', lat, nlon(lat), plon, plev, 1, &
                          qmin(m), arrxzy(:plon,:plev,lat))
            end do

            do k = 1, plev
               call xpavg(arrxzy(:,k,   1), plon)
               call xpavg(arrxzy(:,k,plat), plon)
            enddo

            do k = 1, plev
               do j = 1, plat
                  arrxyz(:,j,k) = arrxzy(:,k,j)
               enddo
            enddo

         end if   ! masterproc

         call scatter_q_field_to_block(arrxyz, m)

      end do
!
!-----------------------------------------------------------------------
! Copy temporary arrays to model arrays
!-----------------------------------------------------------------------
!
      call copy_inidat
!
!-----------------------------------------------------------------------
! Deallocate memory for temporary arrays
!-----------------------------------------------------------------------
!
      deallocate ( ps_tmp )
      deallocate ( arrxyz )
      deallocate ( arrxzy )
      deallocate ( uv_local )
      deallocate ( phis_tmp )        
      deallocate ( landm_tmp )
      deallocate ( sgh_tmp )
      deallocate ( tsice_tmp )
      deallocate ( tsice_rad_tmp )
      deallocate ( tbot_tmp )
      deallocate ( tssub_tmp )
      deallocate ( sicthk_tmp )
      deallocate ( snowhice_tmp )
      deallocate ( pblht_tmp )
      deallocate ( tpert_tmp )
      deallocate ( qpert_tmp )
      deallocate ( landfrac_tmp )
#if ( defined COUP_SOM )
      deallocate ( icefrac_tmp )
      deallocate ( tsocn_tmp )
#endif
!
      return
!EOC
   end subroutine read_inidat
!-----------------------------------------------------------------------


!-----------------------------------------------------------------------
!BOP
! !IROUTINE: copy_inidat --- Copy temporary arrays to model arrays 
!
! !INTERFACE: 

   subroutine copy_inidat 1,24

! !USES:
      use prognostics
      use buffer
      use phys_grid
      use phys_buffer, only: pbuf, pbuf_times, pbuf_get_fld_idx
#if ( defined SPMD )
      use mpishorthand
#endif

      implicit none
!------------------------------Commons----------------------------------
#include <comqfl.h>

! !DESCRIPTION:
!
! Copy temporary arrays to model arrays 
! note that the use statements below contain the definitions
! of the model arrays
!
! !REVISION HISTORY:
!
!   00.06.01   Grant      First attempt at modifying for LRDC
!   00.10.01   Lin        Various revisions
!   00.12.02   Sawyer     Use PILGRIM to scatter data sets
!   01.03.26   Sawyer     Added ProTeX documentation
!   03.08.31   Mirin      Eliminated many 3D temporary global arrays
!
!EOP
!-----------------------------------------------------------------------
!BOC
!
! !LOCAL VARIABLES:
      real(r8), allocatable :: tmpchunk3d(:,:,:)             
      real(r8), allocatable :: tmpchunk(:,:)             
      integer i, j, ic, k, m
      integer n
      integer ncol
      real(r8) :: pmx, pmn

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

      allocate ( tmpchunk(pcols,begchunk:endchunk) )
      allocate ( tmpchunk3d(pcols,plevmx,begchunk:endchunk) )
      tmpchunk(:,:) = 0.
      tmpchunk3d(:,:,:) = 0.
! physics variables
      call scatter_field_to_chunk(1,1,1,plond,landfrac_tmp,landfrac(1,begchunk))
      call scatter_field_to_chunk(1,1,1,plond,landm_tmp,landm(1,begchunk))
      call scatter_field_to_chunk(1,1,1,plond,sgh_tmp,sgh(1,begchunk))
      call scatter_field_to_chunk(1,1,1,plond,tsice_tmp,tsice(1,begchunk))
#if ( defined COUP_SOM )
      call scatter_field_to_chunk(1,1,1,plond,sicthk_tmp,sicthk(1,begchunk))
#endif
      call scatter_field_to_chunk(1,1,1,plond,snowhice_tmp,snowhice(1,begchunk))

      call scatter_field_to_chunk(1,plevmx,1,plond,tssub_tmp,tmpchunk3d)
      do i =begchunk,endchunk
         ncol = get_ncols_p(i)
         surface_state2d(i)%tssub(:ncol,:) = tmpchunk3d(:ncol,:,i)
      end do

#if ( defined COUP_SOM )
      call scatter_field_to_chunk(1,1,1,plond,sicthk_tmp,sicthk(1,begchunk))
      call scatter_field_to_chunk(1,1,1,plond,icefrac_tmp,icefrac(1,begchunk))
      call scatter_field_to_chunk(1,1,1,plond,tsocn_tmp,tsocn(1,begchunk))

! define an initial ocean fraction and non-land ice fraction
! The 1st "where" stmt used to be done in update_srf_fractions (dev45)

      do i =begchunk,endchunk
         ncol = get_ncols_p(i)
         where (icefrac(:ncol,i) + landfrac(:ncol,i) > 1.0)
            icefrac(:ncol,i) = 1. - landfrac(:ncol,i)
         end where

         where (landfrac(:ncol,i) < 1.)
            aice(:ncol,i) = icefrac(:ncol,i)/(1. - landfrac(:ncol,i))
         elsewhere
            aice(:ncol,i) = 0.
         end where
         ocnfrac(:ncol,i) = 1. - landfrac(:ncol,i) - icefrac(:ncol,i)
      enddo
      write(6,*)'INIDAT: ocnfrac=',ocnfrac(1,begchunk)
!
! Master needs global landfrac
!
      call gather_chunk_to_field(1,1,1,plon,landfrac,landfrac_field)
!      write(6,*)'INIDAT iam=',iam,' landfrac=',landfrac
!      write(6,*)'INIDAT iam=',iam,' landfrac_field=',landfrac_field
!
!JR Could read in Focn from initial dataset if available
      Focn(:,:) = 0.
#else
      Focn(:,:) = inf
      frzmlt(:,:) = 0.  ! needs to be 0, otherwise test in tstm always true
      tsocn(:,:) = inf
#endif

! cloud and cloud water initialization should be done in their own packages.  Do it
! here for now since moving it will change answers.

      if (masterproc) then
!
! Arbitrarily initialize all "extra" fields that couldn't be found on the IC file
!
         if(.not. read_pblh) then
            pblht_tmp  (:plon,:) = 0.
            write(6,*) 'Warning:  PBLH not found on IC file; initialized to 0.'
         end if
         if(.not. read_tpert) then
            tpert_tmp  (:plon,:) = 0.
            write(6,*) 'Warning:  TPERT not found on IC file; initialized to 0.'
         end if
         if(.not. read_qpert) then
            qpert_tmp  (:plon,:) = 0.
            write(6,*) 'Warning:  QPERT not found on IC file; initialized to 0.'
         end if
         if(.not. read_tbot) then
!  tbot_tmp now initialized in read_inidat
            write(6,*) 'Warning:  TBOT not found on IC file; initialized with lowest level of T'
         end if
         if(.not. read_tsicerad) then
            tsice_rad_tmp(:plon,:) = tsice_tmp(:plon,:)
            write(6,*) 'Warning:  TSICERAD not found on IC file; initialized with TSICE'
         end if
      endif

      call scatter_field_to_chunk(1,1,1,plond,tbot_tmp,tmpchunk)
      do i =begchunk,endchunk
         ncol = get_ncols_p(i)
         surface_state2d(i)%tbot(:ncol) = tmpchunk(:ncol,i)
      end do
      call scatter_field_to_chunk(1,          1,1,plond,tsice_rad_tmp,tsice_rad(1,begchunk))
      call scatter_field_to_chunk(1,          1,1,plond,pblht_tmp,pblht(1  ,begchunk  ))
      call scatter_field_to_chunk(1,          1,1,plond,tpert_tmp,tpert(1  ,begchunk  ))

! Qpert for only the first constituent is initialized and scattered; for all other
!   constituents, qpert is set to zero.
      call scatter_field_to_chunk(1,          1,1,plond,qpert_tmp,tmpchunk(1,begchunk ))
      qpert(:,:,:) = 0.
      do i =begchunk,endchunk
        qpert(:,1,i) = tmpchunk(:,i)
      enddo
!
! Global integerals
!
      if (masterproc) then
         zgsint = zgsint_tmp
      endif
!
#if ( defined SPMD )
      call mpibcast (zgsint,1,mpir8,0,mpicom)
#endif
      deallocate ( tmpchunk )
      deallocate ( tmpchunk3d)
!EOC
   end subroutine copy_inidat

!-----------------------------------------------------------------------
!BOP
! !IROUTINE: scatter_q_field_to_block --- scatter a 3D constituent array to prognostic array q3
!
! !INTERFACE: 

   subroutine scatter_q_field_to_block(xyz, cnst_idx) 1,3

! !USES:
      use shr_kind_mod, only: r8 => shr_kind_r8
      use pmgrid,       only: plon, plat, plev, beglat, endlat, &
                              beglev, endlev, strip3dxyz
      use prognostics,  only: q3
#if ( defined SPMD )
      use mpishorthand, only: mpicom
#endif

      implicit none

! !INPUT PARAMETERS:
      real(r8), dimension(plon,plat,plev), intent(in) :: &
         xyz        ! 3D constituent field
      integer, intent(in) :: &
         cnst_idx   ! constituent index in prognostic array q3

! !DESCRIPTION:
!
! Scatter a 3D constituent array from the master processor to the 
! q3 array in the prognostics module.
!
! !REVISION HISTORY:
!
!   02.07.31   Eaton      Initial version
!
!EOP
!-----------------------------------------------------------------------
!BOC
!
! !LOCAL VARIABLES:
      integer j, k
      real(r8), allocatable :: xyz_local(:,:,:)

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

#if ( defined SPMD )

      allocate( xyz_local(plon,beglat:endlat,beglev:endlev) )
      call scatter( xyz, strip3dxyz, xyz_local, mpicom )
      do k=beglev,endlev
         do j=beglat,endlat
            q3(:,j,k,cnst_idx) = xyz_local(:,j,k)
         enddo
      enddo
      deallocate( xyz_local )

#else

      do j=beglat,endlat
         do k=beglev,endlev
            q3(:,j,k,cnst_idx) = xyz(:,j,k)
         enddo
      enddo

#endif

!EOC
   end subroutine scatter_q_field_to_block

!-----------------------------------------------------------------------
end module inidat